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1.  Bitroduction 


The  recovery  of  mean  gravity  anomalies,  or  alternate  representations 
such  as  surface  densities,  point  masses,  etc.,  from  the  range  rate  sum  observa- 
tions in  a 'hi^-low'  satellite  experiment  has  been  the  subject  of  several  investiga- 
tions, among  others  in  (Schwarz,  1970;Kaula,  1972;  Martin,  1972;  Hajela,  1974; 
Kaula  et  al. , 1975;  Sjogren  et  al. , 1976;  Vonbun  et  al. , 1977).  The  recovery  of 
gravity  parameters  in  these  investigations  was  based  on  the 'deterministic'  approach 
of  least  squares  approximation,  where  we  do  not  fully  utilize  the  statistical 
information  relating  to  the  parameters,  particularly  the  covariance  functions 
between  the  parameters  and  the  observations  (or  more  usually  a related  quantity 
to  the  observations). 

The  'least  squares  collocation'  approach  (for  a conceptual  discussion, 
see  Moritz,  1976,  Sec.  4)  for  the  recovery  of  gravity  anomalies  from  range  rate 
sum  observations  in  high-low  satellite  to  satellite  tracking  (SST)  was  suggested 
by  Rummel  (1975,  Sec.  5)  utilizing  the  cross-covariance  function  of  gravity 
anomalies  with  the  radial  derivative  (3T/Sr)  of  the  anomalous  potential  (T),  and 
the  autocovariance  function  of  9T/9r.  The  numerical  computation  of  these 
functions  was  made  possible  by  the  covariance  expressions  developed  by 
Tscheming  (1976)  using  an  assumed  model  for  the  anomaly  degree  variances 
(Tscheming  and  Rapp,  1974). 

Rummel,  Hajela  and  Rapp  (1976)  carried  out  simulation  studies  for 
developing  the  procedure  for  the  recovery  of  mean  gravity  anomalies  from  SST 
range-rate  data  using  least  squares  collocation.  The  present  report  applies  the 
results  of  the  previous  study  (ibid.,  1976)  to  the  systematic  investigation  of 
problems  in  utilizing  real  ATS-6/GEOS-3  observations.  These  problems  will  be 
discussed  in  outline  in  Section  1.2  and  in  greater  detail  in  Sections  2 and  3.  The 
recovery  of  eight  5°  mean  gravity  anomalies  in  a local  area  will  be  presented  in 
Section  4. 


1. 1 Recovery  of  Anomalies  from  Simulated  Data 

We  briefly  recapitulate  the  discussion  in  Rummel  et  al.  (1976).  We 
consider  the  range  rate  data  to  be  the  residual  value  referred  to  the  earth's 
normal  gravitational  field  expressed  by  potential  coefficients  (P.C.)  up  to  degree 
and  order  12.  That  is,  we  subtract  from  the  observed  range  rate  sum  (Martin, 
1972;  Hajela,  1974)  at  any  given  time,  usually  at  10  seconds  data  interval,  the 
computed  value  of  the  range  rate  sum  corresponding  to  the  cartesian  position 
(X,  Y,  Z)  and  velocity  (X,  Y,  4)  coordinates  of  ATS-6  and  GEOS-3  satellites. 
These  coordinates,  or  state  vector  (X,  Y,  Z,  X,  Y,  Z),  are  obtained  from  an 
initial  state  vector  at  an  initial  epoch  (usually  a few  hours  earlier)  by  integrating 
the  equations  of  motions  of  the  satellites  in  an  inertial  coordinate  system 
(usually  referenced  to  0**  of  the  day)  in  a force  field  expressed  by  the  earth's 
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low-degree  nomial  gravitation  (12,  12  P.C.).  and  also  taking  into  account  the 
lunar  and  solar  gravitation  and  the  solar  radiation  pressure.  (The  atmospheric 
drag  is  negligible  as  the  altitude  of  satellites  is  higher  than  800  km.)  We  will 
consider  the  determination  of  initial  state  vectors  later,  but  it  is  obvious  that 
the  residual  values  (observed-computed)  of  the  range  rate  would  be  biased  if 
the  initial  state  vectors  are  not  known  accurately. 

The  residual  range  rate  then  expresses  the  range  rate  of  GEOS-3, 
sensed  at  ATS-6,  due  to  anomalous  potential  T 

(1.1)  T = W - U 

where  W expresses  the  actual  gravitational  potential  of  the  earth,  and  U the 
low-degree  (12,  12  P. C.)  normal  potential,  both  at  the  location  of  GEOS-3 
satellite.  The  line  of  sight  (ATS-6/GEOS-3)  residual  range  rate  yields  the  line 
of  sight  residual  acceleration  by  numerical  differentiation.  The  effect  of  any 
errors  in  the  raw  residual  range  rate  data  will  be  discussed  later.  However, 
as  the  residual  accelerations  represent  the  slope  of  the  residual  range  rate 
(function),  they  would  be  much  in  error  if  raw  'noisy'  data  is  used  directly. 

The  projection  of  this  residual  acceleration  on  the  radial  direction 
at  the  GEOS-3  location  (denoted  by  the  subscript  c)  gives  the  radial  derivative 
(ST/9  r)(;  of  the  anomalous  potential,  which  is  used  for  the  recovery  of  residual 
gravity  anomalies  (Ag')  referred  to  the  low-degree  (12,  12  P, C.)  normal 
potential.  We  will  henceforward  generally  assume  the  word  residual  to  be 
implied,  unless  required  for  clarity.  We  will  thus  refer  to  the  residual  line  of 
sight  range  rate  simply  as  range  rate  or  R data.  Analogously,  we  will  generally 
use  the  word  'accelerations'  or  R for  the  residual  line  of  sight  accelerations; 
and  assume  3T/3r  without  the  subscript  c to  refer  to  the  GEOS-3  location. 

We  may  refer  to  residual  mean  gravity  anomalies  simply  as  anomalies.  If  we 
denote  the  angle  at  GEOS-3  at  any  given  time  between  the  radial  direction  and 
the  line  of  sight  ATS-6/GEOS-3  by  (9,  then 


(1.2)  ST/Sr  = R/cos  8 


and  the  predicted  or  recovered  anomaly  Ag  is  obtained  using  least  squares 
collocation  from  a set  of  n 3T/B  r values  by: 

(1.3)  Afe'  (Q)  = C'(Q)^'Tr  (STrTr+  D)"'  (^T/dr) 
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where  Cirtr  is  the  n x n matrix  of  autocovariarwes  <1  n x 1 ~?T /^r  vector  with 
the  noise  matrix  D,  usually  taken  to  be  a diagonal  matrix  uith  the  same  value  for 
each  element  on  the  diagonal:  C (Ql^/tr  is  iht  transixjsed  vector  of  cross-covar- 
iances of  predicted  anomaly  with  the  ?T/9r  values;  and  Q emphasizes  the  fact 
that  the  location  of  the  anomaly  is  on  the  surfac  of  the  earth,  while  the  3T/3r 
values  refer  to  the  GEOS-3  satellite  lcx;ations.  The  |xjint  auto  and  cross-covariance 
elements  are  computed  from  global  covariance  expressions,  and  the  mean  anomaly 
covariance  elements  are  computed  by  numerical  integration  of  the  appropriate 
point  covariances. 

The  standard  deviation  bi^'  of  the  predicted  anomaly  is  computed  by: 


(1.4) 


C^^g'Ir  (Crr  Tr  + D)  ^ C^s' 


where  Co  is  the  variance  of  residual  mean  5°  anomalies  and  other  notations  are  as  in 
equation  (1. 3). 

In  view  of  the  local  nature  of  the  covariance  functions  in  equations 
(1.3)  and  (1.4),  when  referred  to  a degree  12  reference  field  (Rummel  et  al. , 1976, 
p.  15,  Table  2.2),  it  was  found  (ibid..  Sec.  6)  that  we  need  to  consider  3T/3r 
values  only  with  sub-satellite  points  within  7.°5  of  the  center  of  the  5°  anomaly 
blocks.  It  was  also  found  that  it  is  adequate  to  consider  9T/3r  values  at  30 
seconds  data  interval.  These  two  considerations  will  be  retained  in  all  tests  with 
real  data.  The  recovery  of  anomalies  with  real  data  will  be  judged  by  statistics 
comparing  the  magnitudes  of  the  predicted  (£^)  against  the  expected  value,  E (Ag), 
and  by  anomaly  discrepancy  e (Ag),  (see  ibid. , pp.  20  and  31): 


(1.5)  e(Ag')  = - E(Ag') 

(1.6)  E(Ag')  = AgT  - Agpc 

where  Agr  is  the  terrestrial  5°  anomaly  (Rapp,  1977)  and  Agpc  is  the  anomaly 
implied  by  the  degree  12  reference  field. 

The  elements  of  the  diagonal  matrix  ^ in  equations  (1.3)  and  (1.4) 
could  conceivably  be  considered  as  the  square  of  the  standard  deviation  of  9T/^r 
data,  which  may  be  obtained  by  propagation  of  the  variance  of  R , which  in  turn 
could  be  obtained  in  some  manner  from  the  assumed  standard  deviation  of  raw 
R data.  However,  with  both  simulated  and  real  data,  tests  on  the  stability  of  the 
solution  were  performed  with  different  assumed  standard  deviation  of  9T/9r 
data  like  0.5,  1.0,  1.5,  2.0  mgals.  Equation  (1.4)  is  not  very  sensitive  to  this 
variation,  but  in  equation  (1.3),  a standard  deviation  of  2 mgals  or  larger  would 
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show  a dampened  recovery  with  low  R.M.S.  (root  mean  square)  value  of 
predicted  On  the  other  hand,  a small  value  less  than  0.  5 mgals  would 
show  an  unstable  recovery  with  large  RMS  ^ (usually  larger  than  RMS 
value  of  expected  anomalies  E (Ag”^)),  and  also  large  RMS  anomaly  discrep- 
ancy e(Ag).  For  details,  see  Rummel  et  al,  (1976,  pp.  34  and  43). 

Other  details  of  the  formulations  and  procedures  may  also  be 
seen  in  the  above  report,  and  are  not  being  repeated  here.  We  will  now 
describe  the  variations  which  were  required  to  implement  the  use  of  re? 

R data. 


1.2  Recovery  of  Anomalies  from  Real  Data 

In  simulation  studies,  we  considered  the  initial  state  vectors  of 
the  satellites  to  have  been  determined  a-priori  to  an  accuracy  so  that  any 
resulting  error  (modeling  error)  in  the  residual  R data  (due  to  treating  the 
initial  state  vectors  to  be  without  error)  was  much  less  than  the  observational 
errors  or  'noise'.  This  a-priori  knowledge  was  necessary,  as  it  was  found 
that  the  initial  state  vectors  could  not  be  converged  (GEODYN,  Vol.  Ill,  pp.  1. 2 - 
5;  also  see  below)  from  the  limited  time  span  of  30  to  60  minutes  of  range 
rate  sum  observations. 

We  presume  here  that  at  the  satellite  altitudes  exceeding  800  km, 
the  equations  of  motion  of  the  satellites  can  be  expressed  with  negligible  error 
if  the  earth's  gravitational  field  is  described  by  the  full  set  of  potential  coeff- 
icients of  GEM-7  (Wagner  et  al. , 1976).  If  we  also  take  into  consideration  the 
lunar  and  solar  gravitation,  and  account  for  the  ionospheric  and  tropospheric 
refraction  and  the  transponder  delays  at  the  satellites,  then  barring  any  system 
bias  in  the  range  rate  sum  observations,  the  computed  values  should  fit  the 
observed  values  in  the  least  squares  sense,  and  we  may  iteratively  solve  for 
the  initial  state  vectors  of  the  ATS-6  and  GEOS-3  satellites.  We  will  consider 
the  solution  to  have  converged  when  the  corrections  to  the  initial  state  vectors 
in  two  consecutive  iterations  do  not  change  by  more  than  2%. 

The  converged  initial  state  vectors  are  thus  based  on  the  observation 
type  and  the  time  span.  They  may  be  strongly  biased  if  the  observations  were 
biased  in  some  way,  and  this  possibility  would  increase  if  only  one  type  of 
observation  was  used  for  converging  the  arc,  and  also  if  the  time  span  of 
observations  was  'short'.  The  adequacy  of  time  span  depends  on  the  type  of 
observations  (and  their  variance).  Conceptually,  range  computations  will  be 
more  sensitive  to  the  initial  state  vector  of  a satellite  as  compared  to  the  range 
rate  values,  and  the  range  rate  sum  values  of  a high-low  satellite  pair  will  be 
still  less  sensitive  to  the  initial  state  vectors  of  the  satellites.  We  should  there- 
fore expect  the  changes  in  the  solution  of  initial  state  vectors  in  successive 
iterations  to  vary  sharply  with  the  last  observation  type,  and  several  iterations 
may  be  needed  for  the  converged  solution  even  if  the  time  span  was  comparatively 
long,  like  1 to  3 hours. 
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The  optimum  solution  would  then  be  to  use  several  types  of 
observations  obtaining  a weighted  least  squares  solution  of  the  initial  state 
vectors  of  ATS-6  and  GEOS-3.  As  we  did  not  have  access  to  other  types  of 
observations  besides  the  ATS-6/GEOS-3  SSE  tracking  data  (NASA,  1976), 
the  initial  state  vectors  were  not  solved  for  by  us.  These  were  kindly  supplied 
by  Marsh  (1976,  1977)  and  will  be  discussed  in  Section  2.  The  effect  of 
uncertainty  in  the  initial  state  vectors  on  the  residual  R,  R and  anomalies 
will  be  discussed  in  Section  4. 

In  simulation  studies,  the  residual  R data  was  numerically 
differentiated  by  approximating  it  by  a smooth  continuous  function,  viz.  the 
interpolating  natural  cubic  spline,  and  then  analytically  differentiating  the 
spline  to  obtain  R.  The  spline  being  a set  of  piecewise  continuous  cubic 
polynomials  does  not  exhibit  spurious  oscillations  between  the  data  points, 
which  is  a characteristic  of  higher  order  polynomials.  Also,  as  the  spline 
has  the  minimum  norm  among  all  interpolating  functions  (details  in  Section  2), 
it  is  the  smoothest  function  and  its  first  derivative  (R  in  our  case)  is  also  a 
smooth  (and  continuous)  function.  The  interpolating  spline,  however,  passes 
through  the  R data  points  and  does  not  admit  any  errors  in  them.  The  slopes 
of  the  interpolating  spline  would  thus  be  in  much  error  with  real  R data,  which 
have  observational  errors,  as  the  spline  would  have  spurious  ’rise  and  fall’  in 
trying  to  fit  each  data  point  exactly.  We  have  to  thus  consider  a cubic  spline 
which  may  fit  the  data  in  the  least  squares  sense  (Lawson and  Hanson,  1974, 
Chap.  27)  so  that  the  fitted  spline  filters  the  raw  R data  of  the  ’noise’.  A 
similar  data  fitting  approach  was  reported  by  Muller  and  Sjogren  (1968)  in 
filtering  lunar  doppler  data  by  fitting  piecewise  continuous  cubic  Hermite 
polynomials,  but  which  did  not  have  a continuous  f'rst  derivative.  Cubic 
splines  were  used  in  their  later  analyses. 

Besides  filtering  the  raw  R data,  we  also  wish  to  smooth  it  in  the 
sense  that  its  slope,  R,  changes  smoothly  in  a physically  meaningful  way.  We 
know  from  simulation  studies  (Rummel  et  al.,  1976,  p.  54)  that  the  values  of 
ST/9r  (and  R)  change  only  gradually  within  a mgal  or  so  over  30  seconds 
interval  at  the  GEOS-3  altitude.  The  spline  fitted  in  the  least  squares  sense 
should  then  be  further  required  to  be  smooth  enough  so  that  its  first  derivative 
(and  3T/9r)  should  change  in  the  gradual  manner  indicated.  The  smoothing  is 
dependent  on  the  spacing  of  spline  nodes,  where  different  cubic  polynomials  in 
neighboring  intervals  meet.  If  the  spline  nodes  are  too  far,  we  will  get  a 
dampened  representation  of  R from  over-smoothii^  of  the  spline;  while  if  the 
nodes  are  too  near,  the  slopes  of  the  fitted  spline  change  sharply  at  the  data 
points  (at  10  seconds  interval)  giving  noisy  R values.  This  will  be  discussed 
in  detail  in  Section  3. 
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1.3  Scope  of  Investigations 


The  ATS-6/GEOS-3  SSE  tracking  data  was  supplied  by  NASA  (1976) 
in  the  ATSR  format  (Bryan  et  al. , 1975).  A portion  of  this  data  from  April  16, 
1975  to  May  22,  1975  was  reformatted  in  the  GEODYN  binary  format  (GEODYN, 
Vol.  ni,  Appx.  C 7)  using  a preprocessor  program  kindly  supplied  by  Martin 
(1975).  Agajelu  (1977)  selected  five  descending  GEOS-3  revolutions  154,  268, 

439,  453  and  467,  in  which  the  GEC)S-3  satellite  moved  from  north-east  to 
south-west,  for  examining  the  recovery  of  5°  anomalies  in  the  Carribean 
portion  of  the  North  Atlantic  Ocean  between  the  latitudes  15°  to  35°  North  and 
longitudes  275°  to  295°  East.  The  initial  state  vectors  for  both  ATS-6  and 
GEOS-3  satellites  were  supplied  by  Marsh  (1976)  for  these  five  GEOS-3 
revolutions  during  the  period  April  20  to  May  12,  1975. 

It  was  intended  to  use  five  ascending  GEOS-3  revolutions,  in  which 
the  GEOS-3  satellite  moved  over  the  area  of  investigation  from  south-east  to 
north-west.  These  revolutions  were  104,  118,  175,  232  and  246  during  the 
period  April  17  to  April  27,  1975.  The  initial  state  vectors  were  received 
(Marsh,  1977)  only  for  the  GEOS-3  satellite  for  three  revolutions.  It  was 
hoped  that  the  corresponding  initial  state  vectors  for  ATS-6  could  be  obtained 
by  integration  from  the  values  available  for  O*'  on  April  19  and  April  25,  1975. 

The  present  study  investigates  the  above  limited  real  range-rate 
sum  observations  for  establishing  the  procedures  for  the  recovery  of  5° 
anomalies  using  least  squares  collocation.  If  the  range  rate  sum  observations 
in  all  10  GEOS-3  revolutions  could  be  used,  we  would  have  observations  from 
one  ascending  and  one  descending  revolution  for  the  recovery  of  each  5°  anomaly. 
As  we  will  later  see  in  Sections  2 and  4,  this  density  of  observations  was,  how- 
ever, not  available. 

We  will  examine  the  reliability  and  compatibility  of  the  initial  state 
vectors  in  Section  2,  and  their  effect  on  the  recovery  of  anomalies  wQl  be 
investigated  in  Section  4.  We  will  consider  the  filtering  and  smoothing  of  the 
residual  R data  in  Section  3. 


2.  Initial  State  Vectors  for  ATS-6  and  GEC)S-3  Satellites 

The  initial  state  vectors  were  provided  by  Marsh  (1976)  for  both 
satellites  corresponding  to  five  descending  revolutions  of  GEOS-S.  We  first 
describe  in  Section  2. 1 the  particulars  of  the  observations  and  time  span  used 
for  determining  these  initial  state  vectors,  and  the  uncertainties  ascribed  to 
them.  Initial  state  vectors  for  three  ascending  revolutions  of  GEOS-3  were 
later  provided  by  Marsh  (1977),  together  with  data  from  which  it  was  hoped  that 
the  initial  state  vectors  for  two  other  ascending  revolutions  of  GEOS-3,  and  the 
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corresponding  elements  for  ATS-6  (for  all  five  ascending  revolutions  of  GEOS-3), 
may  be  obtained  by  forward  or  backward  integration.  We  discuss  in  Section  2.  3 
attempts  to  recover  these  values.  We  also  describe  in  Section  2.2  the  result 
of  trying  to  filter  out  some  observational  and  modeling  errors  in  a short  time 
span  of  range-rate  observations  by  letting  the  initial  state  vectors  be  constrained 
with  suitable  weights  to  change  only  slightly  to  see  if  some  errors  may  be  absorbed 
by  them.  The  results  showed  that  this  approach  is  not  workable. 


2. 1 Initial  State  Vectors  Available  for  Both  Satellites 


The  five  descending  revolutions  of  GEOS-3  were  154,  268,  439,  453 
and  467.  The  epoch  of  the  initial  state  vectors,  the  time  span  of  observations 
and  other  particulars  are  given  in  Table  2. 1.  The  earth's  gravitational  field 
was  modeled  by  the  full  set  of  potential  coefficients  of  GEM-7  (Wagner  et  al. , 
1976).  The  lunar  and  solar  gravitational  effects,  including  earth  tides,  and  the 
solar  radiation  pressure  were  taken  into  account.  The  atmospheric  drag  was 
considered  zero  at  the  satellite  altitude  exceeding  800  km.  The  observations 
considered  were  range,  including  ’user  ranges,  range-rate  and  average  range- 
rate  (GEODYN,  Vol.  Ill,  1975,  Appx.  C.7).  The  initial  state  vectors  of  both 
ATS-6  and  GEOS-3  were  adjusted  together,  keeping  observations  for  a single 
pass  of  GEOS-3,  as  seen  from  ATS-6,  centrally  in  the  time  span  of  observa- 
tions. The  adjustment  was  done  iteratively  (GEODYN,  Vol.  I,  1976,  Sec.  2.2) 
minimizing  the  sum  of  squares  of  residuals  of  observed  minus  computed  (O  - C) 
values  in  the  force  model,  till  the  RMS  (root  mean  square)  value  in  the  last 
iteration  converged  within  2%  of  the  RMS  value  in  the  previous  iteration.  The 
RMS  position  (RMS  Pi)  standard  deviation  (S  D)  and  RMS  velocity  (RMS  Vi) 
standard  deviation  for  each  satellite  (i  = 1, 2)  was  computed  from  the  standard 
deviations  cf  the  position  (Sxt,  Sn,  Szi)  and  velocity  (Sli,  Sy*i,  SJi)  elements 
at  the  initial  epoch  as: 


(2.1) 

RMS  Pi  = (Sxi^ 

+ Syi^  + 

i = 1 for  ATS-6, 

i = 2 for  GEOS-3 

(2.2) 

RMSVi  = (Sii^ 

+ Syi*'  + 

S?i^)^, 

i = 1 for  ATS-6, 

i = 2 for  GEOS-3 

The  data  interval  for  ange-rate  observations  in  all  cases  was  10 
seconds,  except  for  revolution  467  where  it  was  1 second  during  a period  of 
about  7 minutes  and  10  seconds  otherwise.  There  was  also  a break  of  about 
30  seconds  at  each  end  of  1 second  interval  observations  in  revolution  467. 
The  residuals  in  this  revolution  after  the  convergence  of  initial  state  vector 
were  much  larger  by  a factor  of  4 to  6 compared  to  other  revolutions.  This 
revolution  was  thus  not  used,  and  only  four  revolutions  154  , 268  , 439  and  453 
were  used  in  subsequent  work. 
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Table  2. 1.  Data  Used  for  Determining  Initial  State  Vectors  for  ATS-6  and 
GEOS-3  corresponding  to  GEOS-3  Revolutions  154,  268,  439, 


2.2  Filtering  of  Residual  R Data  by  Slight  Variations  to  Initial  State  Vectors 


FT 


The  statistics  RMSPi  and  RMSVi  , expressing  the  uncertainty  in 
the  determination  of  the  initial  state  vectors,  were  large  (see  Table  2. 1)  for 
all  revolutions,  only  revolution  453  having  a comparatively  lower  value.  This 
was  a cause  for  concern,  for  if  the  initial  state  vectors  were  incorrectly 
determined,  it  will  cause  large  residuals  of  range  rate  (Agajelu,  1976,  Sec.  4). 
A larger  value  of  RMSP  and  RMSV,  however,  does  not  necessarily  imply  a 
large  discrepancy  from  the  'true'  value  of  the  initial  state  vector,  which  is 
ensured  by  taking  different  types  of  observation  for  an  adequate  time  span  and 
by  the  convergence  criterion  of  the  determination  of  the  initial  state  vector. 

We  will  be  using  the  range-rate  R,  residual  to  a (12,  12)  potential 
coefficients  field,  as  the  raw  data  to  recover  residual  gravity  anomalies.  This 
data  is  needed  for  a limited  time  span  for  the  recovery  of  anomalies  in  a local 
area.  We  cannot  use  this  limited  data  to  improve  the  available  initial  state 
vectors  in  the  sense  of  bringing  them  closer  to  the  'true'  values.  But  it  may 
be  argued  that  we  could  poLSibly  filter  the  raw  R values  of  some  observational 
and  modeling  errors  by  letting  the  initial  state  vectors  take  the  'slack'  by 
letting  them  vary  slightly  within  the  standard  deviation  implied  by  the  statistic 
RMSP  and  RMSV.  This  may  still  be  too  large  for  a case  like  revolution  154 
in  Table  2. 1.  There  is  also  a risk  in  this  approach  that  the  initial  state  vectors 
may  absorb  a part  of  the  signal,  while  filtering  the  observational  and  modeling 
errors. 


Nevertheless,  we  tested  this  approach  by  considering  the  case  of 
revolution  453,  where  the  a-priori  values  of  the  initial  state  vectors  as  avail- 
able in  Table  2. 1 were  assigned  standard  deviations  implied  by  items  6 to  9 
(in  Table  2. 1),  and  then  allowed  to  change  to  fit  a span  of  30  minutes  residual 
ft  values.  Three  iterations  were  tried,  and  the  results  are  shown  in  Table  2.2. 

It  is  obvious  that  the  range-rate  observations  by  themselves  are 
not  very  helpful  in  tying  down  the  initial  state  vector.  And  that,  in  any  case, 
a 30  minute  time  span  of  R data  is  not  long  enough  to  do  so.  A shorter  time 
span  would  be  only  worse.  We  have  to  thus  fix  the  initial  state  vectors, 
obtained  from  a longer  time  span  and  with  different  data  types  in  Section  2. 1, 
by  assigning  them  very  low  standard  deviations  (Hajela,  1974,  Section  6. 1). 

The  filtering  and  smoothing  of  observational  and  modeling  errors  would  have 
to  be  handled  subsequently.  This  will  be  discussed  in  Sections  3. 1 and  3.  2. 
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Table  2.2.  Variation  of  Initial  State  Vectors  to  Filter  Residual  Range-Rate 
Data  for  30  Minutes  in  Revolution  453 


1st 

•Iteration 

— 

2nd 

Iteration 

3rd 

Iteration 

1.  Observations  accepted  for  adjustment 

(a)  No.  of  obsns. 

180 

17 

61 

(b)  Mean  value  (cm/sec) 

0.72 

-50.85 

23.07 

(c)  RMS  value  (cm/sec) 

0.82 

67.78 

29.28 

2.  Change  to  Initial  State  Vectors 

(a)  ATS-6  in  position  (meters) 

0.53 

0.40 

0.38 

(b)  ATS-6  in  velocity  (cm/sec) 

0.06 

0.05 

0.07 

(c)  GEOS-3  in  position  (meters) 

1. 19 

12.15 

4.88 

(d)  GEOS-3  in  velocity  (cm/sec) 

0.02 

1.85 

0.08 

3.  Uncertainties  in  Initial  State  Vectors 

(a)  R^BP  for  ATS-6  (meters) 

17.15 

17.21 

17.19 

(b)  RMSV  for  ATS-6  (cm/sec) 

0.44 

0.48 

0.46 

(c)  RMSP  for  GEOS-3  (meters) 

9.90 

12.18 

11.65 

(d)  RMSV  for  GEOS-3  (cm/sec) 

1.05 

1.32 

1.14 

2.3  Initial  State  Vector  Available  for  only  GEOS-3  Satellite 

The  five  ascending  revolutions  of  GEOS-3  (which  we  wished  to 
consider)  with  epoch  times  during  April  75  were:  104  (April  17,  07"),  118 
(April  18,  07"),  175  (April  22,  08"),  232  (April  26,  08")  and  246  (April  27, 
08").  Initial  state  vectors  for  GEOS-3  for  revolutions  175,  232  and  246  were 
received  from  Marsh  (1977).  The  particulars  for  the  observations  used  for 
determining  the  initial  state  vectors  for  GEOS-3  are  given  in  Table  2.3,  on 
the  same  lines  as  in  Table  2. 1.  As  the  tracking  data  before  April  19,  1975 
were  not  readily  available,  the  initial  state  vectors  for  GEOS-3  were  not 
determined  for  revolutions  104  and  118.  Instead,  the  initial  state  vectors 
for  both  ATS-6  and  GEOS-3  satellites  solved  together  were  supplied  for 
April  19  (O"  0®  0®),  and  later  also  for  April  25  (0"  0*  0").  The  particulars 
of  observations  used  for  determining  these  are  also  given  in  Table  2.  3. 

I 

I 

I 
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Table  2.  3.  Data  Used  for  Determining  Initial  State  Vectors  for  GEOS-3 
Revolutions  175,  232,  246  and  for  ATS-6  and  GEOS-3  for 
Epochs  on  April  19  and  April  25,  1975. 


f 
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We  note  that  the  ATS-6  elements  corresponding  to  the  epoch  of 
GEOS-3  revolutions  were  not  available.  It  was  hoped  that  these  could  be 
obtained  by  integration  from  the  ATS-6  elements  available  for  19  April  75, 
or  preferably  from  the  nearer  epoch  25  April  75.  And  similarly,  the  ATS-6 
and  GEOS-3  elements  for  revolutions  104  and  118  could  be  obtained  by  backward 
integration  from  the  eixx;h  on  19  April  75  using  the  GEODYN  (1975)  program. 
Before  this  was  done,  it  was  necessary'  to  examine  the  compatibility  of  the 
elements  supplied  in  the  last  two  columns  of  Table  2.3,  which  had  been 
rigorously  determined  from  observations  of  different  tyjies  over  5 day  periods. 
This  was  first  done  for  the  ATS-6  satellite.  The  true  of  date  ephemeris  for 
ATS-6  was  generated  by  integrating  forwards  and  backwards  in  the  force 
field  as  described  in  Section  2. 1 for  the  period  19  April  75  to  29  April  75 
sep.*rately  from  the  elements  available  in  the  last  two  columns  of  Table  2.3. 

The  position  (X,  Y,  Z)  and  velocity  (X,  Y,  Z)  coordinates  were  compared  from 
these  two  ephemeris,  and  also  with  the  ATS-6  elements  available  in  Table  2. 1. 
The  comparisons  could  be  made  for  the  epochs  on  19  April  0‘,  20  April  for 
epoch  of  revolution  154,  22  April  for  epoch  of  revolution  175,  25  April  0*’, 

26  April  for  epoch  of  revolution  232,  27  April  for  epoch  of  revolution  246,  and 
28  April  for  epoch  of  revolution  268.  Comparisons  were  also  made  with  ATS-6 
elements  supplied  with  SSE  tracking  data  (NASA,  1976)  for  the  epoch  on 
20  April  0*’.  The  differences  are  shown  in  three  rows  for  each  ejxjch  in 
Table  2.4;  (a)  first  row  — integrated  elements  from  April  19,  0^  epoch  minus 
elements  available  in  Tables  2.1  or  2.3;  (b)  second  row  — integrated  elements 
from  April  25,  0**  epoch  minus  elements  available  in  Tables  2. 1 or  2.  3;  third 
row  — integrated  elements  from  April  19,  0*"  epoch  minus  integrated  elements 
from  April  25,  0^  epoch.  When  there  was  no  data  to  complete  a row,  it  was 
left  blank. 


It  is  clear  from  the  entries  against  row  (c)  in  Table  2.4  that  we  may 
be  able  to  integrate  the  elements  reasonably  correctly  only  within  the  time  span 
of  observations  used  for  determining  the  initial  state  vector.  We  may  note  the 
large  values  of  the  discrepancies  (A-B)  in  rows  Ic  and  5c.  But  what  is  more 
surprising  is  that  they  differ  from  each  other  thou^  the  time  period  for  back- 
ward integration  (from  25  April  to  19  April  75)  of  ATS-6  elements  on  25  April  75 
is  the  same  as  the  time  period  for  forward  integration  (from  19  April  to  25  April  75) 
of  ATS-6  elements  on  19  April  75.  We  may  also  note  the  very  large  changes  in 
rows  7c  and  8c,  and  also  the  changes  in  rows  4c  and  5c. 

The  ATS-6  elements  on  19  April  75  and  25  April  75  were  converged 
independently,  based  on  5-day  observations  with  low  RMSP  and  RMSV  values 
(Table  2.3).  But  they  are  not  at  all  compatible  with  each  other  outside  the  time 
span  of  observations  on  which  they  are  based.  This  is  also  seen  in  the  large 
discrepancy  in  row  2a  of  Table  2.  4.  The  elements  differ  widely  as  they  are 
based  on  different  time  spans  of  observations.  The  differences  are  grossly 
large  in  row  8b.  We  recall  that  the  converged  elements  (C)  for  ATS-6  in  row 
8b  were  obtained  for  time  span  of  about  2 1/2  hours  (see  Table  2. 1 for 
revolution  268). 
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Table  2.4.  Differences  in  Position  and  Velocity  Coordinates  for  ATS-6  as  Integrated 
from  Epochs  on  19  April  75  (A)  and  25  April  75  (B)  in  Comparison  with 
Converged  Elements  (C)  for  Other  Epochs 
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It  is  therefore  not  possible  to  integrate  backwards  from  the  initial 
state  vectors  of  ATS-6  and  GEOS-3  on  19  April  75  to  obtain  elements  corresponding 
to  GEOS-3  revolutions  104  and  118  on  17  April  and  18  April  75  respectively.  We 
may  also  integrate  the  ATS-6  elements  forward  only  within  the  time  span  of 
observations  used  for  converging  then,  i.  e.  for  epoch  of  revolution  175  from 
elements  on  19  April  75  and  for  epochs  of  revolutions  232  and  246  from  elements 
on  25  April  75.  It  was  found  that  there  was  a break  in  range-rate  observations 
in  GEOS-3  revolution  175  over  the  area  of  investigation,  so  we  could  only  process 
the  observations  in  GEOS-3  revolutions  232  and  246. 

We  then  examined  the  residuals  of  range-rate  sum  in  the  full  force 
field  of  Section  2. 1,  using  the  initial  state  vectors  of  ATS-6  as  obtained  by 
integration  from  25  April  75,  along  with  the  converged  initial  state  vectors  of 
GEOS-3  as  available  from  Table  2.  3,  for  epochs  of  revolution  232  and  246;  56 
and  61  minutes  of  range-rate  sum  observations  were  used  respectively  for  the 
two  revolutions.  We  would  expect  the  residuals  to  be  random  and  close  to  zero, 
if  the  initial  state  vectors  were  satisfectory.  However,  this  was  not  so,  partic- 
ularly for  revolution  232,  and  is  shown  in  Table  2.5.  As  the  GEOS-3  elements 
were  already  converged  as  in  Table  2.  3,  only  ATS-6  integrated  elements  could 
be  responsible  for  unsatisfactory  residuals.  We  then  fixed  the  initial  state 
vector  of  GEOS-3  by  assigning  very  low  standard  deviations,  and  allowed  the 
56  and  61  minutes  of  range-rate  observations  to  iteratively  improve  the  ATS-6 
elements  to  fit  these  observations.  The  a-priori  standard  deviation  of  each 
ATS-6  element  was  taken  approximately  equal  to  the  magnitude  of  the  differences 
(A-B)  in  rows  6c  and  7c  of  Table  2. 4.  A total  of  up  to  9 iterations  were  tried 
separately  for  revolutions  232  and  246,  and  the  results  are  given  in  Table  2.  5. 

We  first  note  that  range-rate  sum  observations,  even  the  time-span 
of  60  minutes,  are  not  very  helpful  in  tying  down  the  initial  state  vector.  (Also 
see  Eddy  and  Sutermeister,  1975,  Sec.  3.)  The  elements  did  not  converge  even 
after  9 iterations  in  revolution  246,  and  even  in  the  case  of  revolution  232  there 
are  very  sharp  fluctuations  in  the  second  and  fourth  iterations.  Secondly,  the 
large  initial  values,  shown  under  iteration  1,  of  mean  and  RMS  residuals  in 
revolution  232,  as  compared  to  revolution  246,  points  to  some  systematic  errors 
in  range-rate  sum  observations  in  revolution  232.  Hence,  even  if  the  elements 
converge  after  seven  iterations,  and  reduce  the  magnitude  of  residuals,  the 
systematic  error  may  still  remain.  This  will  be  discussed  again  in  Section  4.2. 

For  future  reference  we  will  refer  to  the  initial  state  vectors  of  GEOS-3  in 
revolution  232  as  in  Table  2.3  and  the  corresponding  initial  state  vector  of  ATS-6 
as  referring  to  revolution  232  I with  integrated  elements  of  ATS-6,  and  as 
revolution  232  A with  adjusted  elements  of  ATS-6  after  7 iterations  in  Table  2.5. 

The  initial  state  vectors  of  ATS-6  for  revolutions  232  I and  232  A are  shown  in 
Table  2.6.  The  difference  between  them  is  as  grossly  large  as  for  the  case  in 
row  8b  in  Table  2.  4. 
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Table  2.5  Residuals  of  Range-Rate  Sum  in  Full  GEM-7  Field  for  Revolutions  232 
and  246,  Using  Integrated  Elements  of  ATS-6  from  Epoch  25  April  75, 
O'* , and  Convergence  of  ATS-6  Elements. 


(b)  Y (km)  -0.67  0.73  0.55  5.1  8.4  -6.7  22.1  -.018  -1.35 

(c)  Z (km)  -0.79  0.99  -2.09  2.9  -3.8  4.3  6.2  -.020  -1.65 

(d) X(cm/sec)  -2.9  2.4  2.8  40.6  30.8  -19.2  157.4  -0.49  -12.7 

(e)  Y (cm/sec)  4.0  -14.9  6.2  -38.1  -57.6  45.7  -161.1  0.23  10.1 

(f)  Z (cm/sec)  -9.9  20.8  -7.2  -99.7  37.7  -59.6  -295.9  0.80  45.9 


Table  2.6.  Initial  State  Vectors  for  ATS-6  for  Revolution  232 
on  26  April  75  07" 


1 

Revolution 

X 

(m) 

Y 

(m) 

Z 

(m) 

X 

(cm/sec) 

Y 

(cm/sec) 

■ 

Z 

(cm/sec) 

232  I 

-30,047,830 

-29,592,060 

-433,970 

2,156.92 

-2,189.72 

45.05 

232  A 

' Difference 

-30,076,114 

-29,573,587 

-424,089 

2,158.12 

-2,190.88 

41.83 

, (A- I) 

28,284 

-18,473 

-9,881 

-1.20 

1.16 

3.22 

The  convergence  of  the  initial  state  vector  by  itself  is  therefore  not 
enough.  We  have  to  consider  the  observation  type  as  well  as  the  time  span  of 
observations  used  for  the  convergence.  The  integrated  elements  of  ATS-6  are 
totally  discordant  from  the  converged  elements  in  revolution  232  A,  as  well  as 
in  revolution  268  (Table  2.4).  We  will  later  find  in  Section  4.  3 that  the  converged 
elements  of  GEOS-3  in  revolution  268  were  also  grossly  different  from  the 
integrated  elements  from  the  epoch  o''  on  25  April  75.  Firstly,  this  may  point 
to  some  systematic  errors  in  the  data  during  the  time  span  of  about  2 1/2  hours 
used  in  Table  2. 1 for  convergence  of  revolution  268.  Secondly,  there  is  a 
possibility  while  converging  the  elements  in  a satellite  pair,  that  two  widely 
differing  sets  of  initial  state  vectors  may  fit  the  same  range-rate  sum  and  range 
sum  observations.  This  may  be  avoided  by  including  sufficient  independent  range 
observations  to  each  satellite. 

We  will  also  find  later  in  Section  4 that  the  residual  R data  in  revolution 
232  and  revolution  268  were  not  usable  whether  we  use  the  integrated  elements  or 
the  converged  elements  in  the  initial  state  vectors  of  the  two  satellites. 

3.  Computation  of  the  Radial  Derivative  of  the  Anomalous  Potential 

We  consider  in  this  report  the  residual  range-rate  (R)  as  the  raw  data, 
which  is  range-rate  sum  observations  in  satellite  to  satellite  tracking  reduced  by 
the  computed  value  corresponding  to  the  gravitational  field  defined  by  (12,  12) 
potential  coefficients  (Hajela,  1974)  with  the  initial  state  vectors  obtained  in 
Section  2.  In  Rummel  et  al.  (1976),  a natural  cubic  interpolating  spline  was  used 
to  fit  simulated  R data,  which  was  free  of  observational  and  modeling  errors. 

As  the  R data  in  this  report  has  observational  noise,  this  has  to  be  first  filtered 
out  before  further  processing.  We  will  discuss  the  use  of  cubic  splines  in  fitting 


noisy  data  in  the  least  squares  sense  in  Section  3. 1.  The  results  of  this  fitting 
(filtering  and  smoothing),  and  its  analytical  differentiation  to  obtain  residual 
acceleration  (R)  will  be  presented  in  Section  3.2.  The  computation  of  the 
radial  derivative  of  the  anomalous  potential  (9T/Sr)  will  be  described  in 
Section  3.  3. 


3.1  Use  of  Cubic  Splines  for  Fitting  Data 

Let  a set  of  n distinct  and  monotonically  increasing  points  be  defined 
on  the  real  line  on  a closed  interval  I = [a,  b],  such  that; 

(3.1)  a^t^^ts'*^...  tji^b 

and  let  a real  valued  continuous  function  R = f (t)  assume  the  values; 

(3.2)  Ru  = f(tK)  at  the  points  tu  , k = 1,  . . . , n 


We  approximate  this  function  f over  I by  a smooth  'spline'  function 
S = S(t)  consisting  of  a set  of  piecewise  polynomials  of  degree  2r  + 1 defined 
respectively  in  the  intervals  I|(  — ti,  < t ^ t|(,.  I , such  that  S is  2r  times 
continuously  differentiable  at  each  of  the  'nodes'  (or  'breakpoints',  or  'joints') 
ti,  . This  spline  function  is  of  degree  2r  + 1.  We  will  consider  in  this  report 
only  a cubic  spline  (function)  with  r = 1,  which  is  twice  continuously  differen- 
tiable at  each  node  and  is  a set  of  piecewise  cubic  polynomials  in  each  interval 
Ik  . If  we  specify  the  value  of  the  spline  S at  the  nodes  tk  to  be  the  same  as 
the  function  f,  which  is  being  approximated  over  I,  i.e. 

(3.3)  Sk  = S(tk)  = f(tk)~  Rk  » k — l,...,n 


then  S is  the  interpolating  cubic  spline;  and  is  uniquely  defined  if  we  also 
specify  two  additional  boundary  conditions  (Ahlberg,  Nilson  and  Walsh,  1967, 
page  11).  The  dimension  of  S|  is  therefore  n + 2,  where  Si  is  the  linear 
space  of  all  cubic  splines  on  I,  with  nodes  at  tk  , k = 1,  . . . , n;  i.e. 

(3.4)  Si  = [q(t)  e (I)  | q (t)  is  a cubic  polynomial  on 

each  subinterval  U = tk  < t s tk  + i , defined 
on  the  interval  I ] 
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1 


I 

1 


I 


and  (I)  denotes  all  twice  continuously  differentiable  functions  on  I.  Si 
obviously  includes  the  natural  interpolating  cubic  spline  used  in  Rummel 
et  al.  (1976),  but  we  are  now  interested  in  other  splines  also  with  the 
common  condition  that  they  all  have  nodes  at  t^  . 

A unique  representation  for  the  cubic  spline,  and  we  may  refer  to 
it  henceforward  simply  as  spline,  in  terms  of  n + 2 linearly  independent 
members  of  Si  is  thus: 


n 3 

(3.5)  S (t)  = Y Cj  Qj  (t)  , j = 1,  . . . , n + 2,  Cj  is  a real  number. 

j = 1 

If  the  sample  values  of  the  function  f at  the  nodes  tu  are  not  known 
accurately  in  equation  (3.2),  say  due  to  observational  errors  in  Rk»  the 
approximation  of  f by  a spline  using  the  interpolating  condition  (3.  3)  will  be 
unsatisfactory.  The  errors  will  be  further  magnified  in  subsequent  operations 
on  the  spline,  for  example,  if  the  differential  operator  is  of  interest,  then 
S (t)  may  not  recover  the  signal  or  the  trend  (f  (t))  from  noisy  data  Rj,  . 


We  may  then  fit  a spline  S*  (t)  to  the  function  f in  the  least  squares 
sense,  such  that  the  euclidean  or  the  L»  norm  Hf  - S*  lU  is  minimized, 
(Schultz,  1973,  Chap.  6;  Rice,  1969,  Vol.  II«  Sec.  10.4),  i. e.: 


(3.6) 


II  f - S* 


) J (f(t)  -S*(t))^  dt 

_J  w tl( 


n — i 

= I J 


tk+ 1 
tk 


(Ri  - S*(t,))^  dt 


minimum 


and  the  integral  may  be  replaced  by  a summation  for  the  available  observations 
Ri  in  the  interval  I],  . We  now  assume  an  available  vector  R of  noisy 
observations; 


(3.7) 


Rj  = f(ti),  i=l,  ...,  m,  m>n+2 


Following  (3.5)  we  may  represent  S*(t)  as: 


(3.8)  SMt)  = y Cj*  qj  (t)  , j = 1,  ...,  n + 2 

J=  1 


I 
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where  qj  are  (n  + 2)  cubic  polynomials  being  the  bases  of  S|  in  equation  (3.4), 
and  c*  is  the  vector  of  (n  + 2)  unknown  coefficients  Cj*. 

We  then  have  to  solve  for  c*  under  the  constraint  (3.  6)  from  a set 
of  m linear  (observation)  equations: 


(3.9)  » where 

(3. 10)  ^ = [a,j]  ; aij  = qj  (tj)  ; i = 1,  . . . , m;  j = 1,  . . . , n + 2 

There  exist  bases  for  Si  with  the  property  that  if  the  data  are 
ordered,  i.  e.  t^^ts^  ...  - t,,,  then  the  (m  x (n  + 2))  matrix  A is 
band  limited  with  a bandwidth  of  4 (Lawson  and  Hanson,  1974,  page  223). 

We  first  note  that  the  nodes  of  the  spline,  i.  e.  tjj  in  equation  (3. 1) 
are  embedded  in  the  ordered  data  ti  with: 


(3. 11)  tit  - 1 1 for  k = 1,  i = 1 and  for  k = n,  i = m; 

m > n + 2 


We  then  define: 

(3.12)  ^fk~tk  + i~tk»k  = l,...,n-l 

and  scale  the  independent  variable  ti  to  Uj  in  each  interval  1^,  such  that 


1 1 - tk 
^k+l  “ tk 


1 - U 1 


tk  + 1 ~ 1 1 _ tk4.i.  - ti 

t k +1  ~ t k At  k 


Then  for  each  interval  Ik  only  4 cubic  polynomials  out  of  qj(t)  have 
non-zero  values,  i.  e.  for 


(3.14)  j=k-l  + £;  k = l,  ...,n-l;  2 = 1 4; 

qj  (t)  = Pi(l-u) 

qj  + i(t)=  P2(l-u) 

qj  + s(t)=  Ps(u) 
qj+3(t)=  Pi(u) 
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where  (ibid. , 1974)  suggests  the  choice: 


(3.15)  Pi (u)  = 0.25  , and 

P2  (u)  = 1 - 0.75  (1  + u)  (l  - u)^ 


The  spline  at  any  point  t j is  then  evaluated  using  equation  (3.  8)  as 

4 

(3.16)  S*(ti)=y  qk-i+^  (ti);  k = 1,  . . . , n - 1 

£=i 


It  is  easy  to  see  from  equations  (3. 13)  to  (3. 15)  that  S*  is  a cubic 
polynomial  in  each  interval  . It  can  be  shown  that  it  is  continuous  at  the 
nodes  also,  for  in  the  interval  I|(  , when: 


(3.17) 


1 1 


tk 


Uk  = 0 , 


1 - Uk  = 1, 


Qk  (tk)  = 
Qk  + l (f  k)  “ 
qk  + 8(tk)  ~ 
qk+  (f  k)  ~ 


Pi  (1)  = 0.25 

P2(l)  = 1 

Pa  (0)  = 0.25 
Pi  (0)  = 0 


and  when 


(3.18) 


1 1 ~ f k+  1 


Uk+i-1,  l-Ujc+l 


0, 


Qk  (tk+l)  = 
qk  + 1 (tk+l)  = 

Q k+  S (tk+  i)  = 

qk+  3(tk+  l)  ~ 


Pi(0)  =0 
P2(0)  = 0.25 
P2(l)  = 1 
Pi(l)  = 0.25 


Then  if  we  use  the  notation  S*^(ti(+i)  for  the  value  of  spline  in  the 
interval  Ik  at  the  node  Ik+i,  and  similarly  use  the  notation  Sk+i  (tk+i)  for 
the  value  of  spline  in  the  interval  Ik+i  at  the  same  node  Ik+i  , then  the 
continuity  of  S*  is  ensured  by  the  condition: 


(3.19) 


Sk  (tk+l)  " S k + i(tk+l) 


which  is  true,  as  from  equations  (3.16)  to  (3.18): 


-20- 


(3.20)  S*  (tk+i)  = c*  • Pi(0)  • pa(0)  + c f+a  ' P2(l)  + cj+g  • pi  (1) 

~ 0a25c’|t-f)^  + ^ic^s  ^ 0*25c  *+  3 
S*+i(tK+i)  = c *+i  • pi(l)  + cif+a  • pa(l)  + c?+3  • Pa  (0)  + ci!‘+4  • pi(0) 
= 0.25c?+i  + cj'+a  + 0.25  c*+3 

Similarly,  we  may  show  that  the  first  and  second  derivatives  of  the 
spline,  S*'  and  S*  are  also  continuous  at  the  nodes.  From  equation  (3.15): 

(3.21)  Pi"  (u)  = 0.  75  u®  ; Pi"(0)  = 0,  pi'(l)  = 0.  75 

Pa"(u)  = 0.75  (1  + 2u  - 3u^);  p2"(0)  = 0.75,  P2'(l)  = 0,  and 

(3.22)  Pi""(u)  = 1.5u  ; Pi""(0)  =0,  Pi""(l)  = 1.5 
Pa""(u)  = 1.  5 (1  - 3u);  pa"  (0)  = 1.  5,  P2"(l)  = -3 

Using  notations  similar  to  that  in  equation  (3. 19)  and  (3.20),  and 
using  equations  (3.21)  and  (3.22),  it  follows  that: 


(3.23)  sTlt  k + l)  ~ Ck  • Pi(0)  + C)(  + i *Pa(0)+C*+2*P2(l)+Ci(+3  'Pi  (1) 

= 0.75  Ck+1  + 0.75  c*+3 

Sk  + i(tk  + i)  = c • Pi'(l)  +0^+3'  P2"(l)  + cif+a  • Ps"  (0)  + Ck+4  • Pi  (0) 
= 0.  75  c k + i + 0.  75  c k + 3 f and 


(3.24)  S *"(tk  + i)  = c*  • Pi""(0)  + Ck+i*  P2""(0)  + c *+a  • pa  (1)  + c*+3  • Pi"  (1) 

— 1.5  Ck+i“5Ck  + 3 1.5Ck+3 

Sk  + l(tk+l)  = Ck  + i*Pi(l)  + Ck+a'P2(l)  + Ck  + 3 • P2(0)  + Ck  + 4 •Pl(0) 
■”la5cj{^^  “ 1#5ci{  + 3 


1 
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which  leads  to: 


(3.25)  S|(  (t|(+i)—  Si(  + i(tic+i)  • and 

Sic  (tlc+l)  = Sk+l(ti(+i) 

Finally,  it  can  be  shown  that  the  polynomials  in  equations  (3. 17) 
and  (3.18)  do  satisfy  the  condition  (Spath,  1974,  page  58): 


(3.26)  det 


Qk  (tk) 
Qk  (f  k+l) 
qk  (tk) 
q/  (tk+i) 


qk+i  (tk) 
qk  + l (tk+l) 
qk  + l (tk) 
qk+i  (tk+i) 


qk  + S (tk  ) 
qk  + S (tk+l) 
qk  + a (tk ) 
qk  + 2 (tk  + l) 


q k+  3(tk ) 
qk  + 3(tk+ 1) 
qk+ 3(tk) 

qk  + 3(tk+i) 


4 0 


k = 1, 


• • • * 

n - 1 


because  using  equations  (3.17),  (3. 18)  and  (3.21),  equation  (3.26)  takes  the 
following  form: 


det 

Pi  (1) 

P3(l) 

Pa  (0) 

Pi  (0) 

Pi  (0) 

Pa  (0) 

Pa  (1) 

Pi  (1) 

Pi  (1) 

Pa  (1) 

Pa  (0) 

pi  (0) 

Pi  (0) 

Pa'  (0) 

Pa  (1) 

pi  (1) 

det 

0.25 

1 

0.25 

0 

0 

0.25 

1 

0.25 

0.75 

0 

0.75 

0 

0 

0.75 

0 

0.75 

Eqi  ation  (3. 16)  using  equations  (3. 12)  to  (3. 15)  is  therefore  a valid 
represent  tion  of  the  cubic  spline  (3.  8).  We  now  return  to  the  solution  of 
(n  + 2)  vector  £*  from  m linear  equations  (3.9)  under  the  constraint  (3.6). 
This  was  done  by  an  algorithm  given  in  Lawson  and  Hanson  (1974,  Chap.  27). 
It  is  known  that  there  exists  a (mx  m)  orthogonal  (Householder's)  matrix 
Q such  that  if  it  multiplies  from  the  left  the  augmented  (m  x (n  + 3))  matrix 
["a  : of  equation  (3.9),  the  augmented  matrix  is  decomposed  into  a 

((li  + 2)  X (n  + 2))  upper  triangular  matrix  R in  the  following  form: 
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(3.27) 


a Q B [ a A n+3  I a ] — 


n+a^n+2  : n+2 

1 ^n+a  ; 1^1 

B-n-a^n+a  : a-n-a 


and  the  upper  triangular  matrix  R retains  the  same  bandwidth  4 as  the  matrix  A. 
Further,  the  solution  of  (n  + 2)  linear  equations:  ~ 


(3.28) 


a + aRn+a  n+a  - n+adi 


is  the  same  as  that  of  equation  (3.9)  under  rhe  constraint  (3.  6),  and  the  euclidean 
or  La  norm  of  the  residual  vector  r is  given  by  H of  equation  (3.27),  i.e. , 


(3.29) 


Hr  lla  = Hr  - A c^  lU  = I 


The  solution  of  c*  in  equation  (3.28)  is  straightforward  by  back 
substitution  as  R is  upper  triangular.  We  therefore  transform  the  observation 
equations  (3.9)  directly  by  the  QR  decomposition  into  equations  (3.28),  and  the 
formation  of  normal  equations  and  the  inversion  of  normals  matrix  is  avoided. 
Further,  as  there  are  only  five  non-zero  elements  in  the  augmented  matrices 
[A  : R]  and  [^:  d]  in  equation  (3.27),  we  only  need  an  array  (m  x 5)  in  the  storage. 
The  Lawson-Hanson  algorithm  (ibid,  1974)  processes  the  observations  sequentially 
in  this  array  according  to  equation  (3.27),  with  all  observations  in  the  interval 
lit  being  processed  at  a time. 

Specifically,  the  matrix  Q is  a product  of  (n  + 3)  orthogonal  (m  x m) 
Householder's  transformation  matrices  Qj : 


(2.  30)  2 3 ‘ Si+i  2^  * * • 2 ® 2 ^ 

which  reduce  the  column  vectors  ^ ; j = 1,  . . . , n + 3 of  the  augmented  matrix 
[A  : R]  to  0 below  the  diagonal.  The  matrix  ^ j are  computed  by: 

(3.31)  gj  = ^a  - 2 Vj  VjV(v/  Vj) 

where  _I,  is  the  (m  x m)  unit  matrix  and  vector  Vj^  is  the  transpose  of  vector 
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3.2  Filtering  and  Smoothing  of  Raw  R Data 


The  spline  in  equation  (3.16)  filters  the  given  raw  data  R in  the  least 
squares  sense.  It  also  smooths  the  raw  data  depending  on  the  spacing  of  nodes  in 
equation  (3. 12),  or  analogously  depending  on  the  number  of  nodes  in  equation  (3. 1). 
For  a given  number  of  data  points  m,  if  the  number  of  nodes  n = m - 2,  we  get 
the  interpolating  spline.  But  as  n is  reduced,  i.e.  the  spacing  of  nodes  is 
increased,  we  get  a smoother  representation  of  the  data.  It  is  important  to 
choose  a suitable  spacing  of  the  nodes,  particularly  if  S*  (t)  is  of  interest,  as 
the  slope  of  the  fitted  spline  is  sensitive  to  the  smoothing  of  the  data.  Inadequate 
smoothing  will  give  large  and  'noisy'  slope  values,  while  over-smoothing  will 
dampen  them  out.  We  now'  describe  these  exp)eriments  with  different  spacing  of 
nodes  in  respect  to  residual  range-rate  ( R)  data.  The  spacing  of  nodes  was  kept 
uniformly  as  30,  40,  60,  80  or  100  seconds. 

We  first  consider  the  range-rate  sum  observations  in  revolutions  154 
and  453  described  in  Section  2. 1.  We  deduct  from  these  the  computed  value  of 
range-rate  sum  in  the  gravitational  field  of  (12,  12)  potential  coefficients  out  of 
GEM-7  (Wagner  etal.,  1976),  with  the  initial  state  vectors  kept  fixed  at  the  values 
given  in  Table  2.1.  We  consider  these  residual  range -rate  R values  as  the  raw 
data.  We  fit  splines  in  the  least  squares  sense  to  this  raw  data  at  10  seconds 
interval  using  equation  (3. 16)  over  a 20  minute  period,  i.  e.  121  raw  data  points. 
The  20  minute  period  w'as  chosen  to  straddle  the  observations  over  the  area  of 
investigations  centrally.  With  121  raw  data  points,  and  with  the  condition  in 
equation  (3. 11),  the  spacing  of  nodes  at  any  of  the  values  30,  40,  60,  80  or  100 
seconds  ensured  that  the  spline  nodes  were  always  coincident  with  the  location  of 
some  data  points. 

We  give  in  Table  3. 1 for  each  nodal  spacing  the  RIMS  values  in  cm /sec 
of  raw  data,  smoothed  data,  residuals  after  spline  fit  (=  smoothed  - raw  data)  and 
also  the  total  number  of  residuals  out  of  121  points,  w'hich  had  an  absolute  value 
larger  than  an  arbitrary’  value  of  0. 1 cm/sec. 

A plot  of  raw  and  smooth  data  corresponding  to  nodal  spacing  of  40, 

60  and  80  seconds  is  shown  in  Figures  3.1  to  3.3  for  revolution  154,  and  in 
Figures  3.4  to  3.  6 for  revolution  453.  From  Table  3. 1,  as  well  as  Figures  3. 1 
to  3,6,  we  note  greater  smoothing  of  data  with  increased  spacing  of  nodes.  We 
note  that  the  outlying  data  points  have  a lesser  effect  on  the  shape  of  the  spline 
as  the  spacing  of  spline  nodes  is  increased.  This  may  be  seen  particularly  for 
a large  positive  data  point  in  the  center  of  Figures  3. 1 to  3.  3,  and  for  a negative 
data  point  to  the  left  in  Figures  3.4  to  3,6.  The  least  squares  nature  of  fit  is  seen 
clearly  in  Figures  3.  7 and  3,  8,  showing  the  plot  of  residuals  after  spline  fit  at 
nodal  spacing  of  60  seconds  for  revolutions  154  and  453  respectively.  The  RMS 
value  of  residuals  (filtering)  does  not  show'  large  changes  in  Table  3. 1,  but  the 
slopes  (smoothing)  of  the  fitted  splines  get  smoother  (Figures  3. 1 to  3.  6). 
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This  is  seen  more  clearly,  when  we  differentiate  the  spline  (function) 
analytically,  and  plot  S*  , which  are  the  residual  acceleration  R.  The  plots  of 
ii  for  revolution  154,  for  spline  nodes  at  30  , 40  , 60  , 80  and  100  seconds,  are 
shown  in  Figures  3.9  to  3. 13.  We  note  the  'noisy'  oscillatory  pattern  in  Figure  3.  9 
w'hen  the  spline  nodes  are  not  far  enough  at  30  seconds  spacing  to  achieve  adequate 
smoothing  of  data.  On  the  other  hand,  when  the  nodes  are  too  far  at  100  seconds 
spacing,  the  accelerations  are  dampened  off  in  Figure  3. 13.  The  spacing  of  nodes 
at  60  seconds  appears  to  be  optimum  to  filter  and  smoothen  the  raw  R data,  when 
we  are  interested  in  ii.  We  cannot  notice  a clear  difference  between  40,  60  and 
80  seconds  spacing  of  nodes,  when  we  look  at  smoothed  R data  in  Figures  3. 1 to 
3.6  or  in  Table  3. 1.  But  the  Ii  values  do  show  an  optimum  recovery  of  signal 
at  60  seconds. 

This  is  seen  again  in  Figures  3. 14  to  3. 18,  which  shown  the  ii  values 
for  spline  node  spacing  at  30  , 40  , 60  , 80  and  100  seconds  for  the  least  squares 
fitting  of  raw  R data  in  revolution  453.  The  spacing  of  nodes  at  30  seconds  and 
100  seconds  are  clearly  unsuitable,  and  spacing  of  60  seconds  appears  to  be 
optimum  in  the  sense  of  ii  values  changing  gradually  as  we  would  expect  them  to 
at  about  850  km  height  (see  3T/5r  values  in  simulation  studies  (Rummel,  et  al. , 
1976,  page  54)). 

We  also  notice  the  spurious  values  of  accelerations  clearly  in  Figures 
3. 11  and  3. 16,  and  also  in  other  figures  showing  plot  of  accelerations,  near  the 
ends  of  the  data  span.  We  do  not  expect  the  fitting  of  the  spline  to  be  very  satis- 
factory at  either  end  of  the  data  span,  and  this  is  seen  more  clearly  in  the  slopes 
of  the  spline,  than  in  the  spline  itself  in  Figures  3. 1 to  3.6.  With  the  data  interval 
as  10  seconds,  the  data  span  of  20  minutes,  straddling  the  area  of  investigation 
centrally,  was  purposely  chosen  to  be  much  larger  than  required.  2 to  4 minutes 
of  accelerations  appear  to  be  spurious  at  either  end  of  the  data  span,  leaving  at 
least  12  minutes  of  useful  data.  6 to  8 minutes  of  this  data  taken  centrally  was 
actually  used  in  the  recovery  of  anomalies  in  Section  4. 

Information  about  the  RMS  and  mean  values  in  cm/sec  of  the  raw  R 
data,  smoothed  R data  after  fitting  a spline  with  spacing  of  nodes  at  60  seconds, 
and  the  residuals  (smoothed  - raw  data)  is  given  in  Table  3.2  for  all  revolutions 
available  in  Section  2. 1. 

The  large  RMS  value  of  the  residuals  in  revolution  467,  almost  as 
large  as  the  raw  data  itself,  shows  that  these  observations  are  unusable  for  further 
processing.  The  location  of  GEOS-3  in  latitude,  longitude  and  height  in  revolutions 
467  and  268  were  matched  quite  closely,  and  we  would  expect  similar  R values. 

The  RMS  value  of  raw  R data  in  revolution  268,  however,  appears  to  be  much 
larger.  The  plots  of  R and  R in  revolution  268  with  a spline  node  spacing  of  60 
seconds  are  shown  in  Figures  3. 19  and  3.  20.  Figure  3. 19  is  on  a reduced  scale 
of  1/3  as  compared  to  Figures  3. 1 to  3.6  to  accommodate  much  larger  R values. 
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Spline  Nodes  30  sec. 
Data  Interval  10  sec 
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Revolution  453  with  Spline  Nodes  30 
sec.  apart.  Data  Interval  10  sec. 


+ 8 mgals 


8 mgals 


+ 8 mgals 


- 8 mgals 


8 mgals 


8 mgals 


8 cm/sec 


cm/sec 


+ 8 mgals 


- 8 mgals 


t 


Table  3.2.  Least  Squares  Spline  Fit  with  Spline  Nodes  60  Seconds 


Apjart  to  Raw  R Data  at  10  Seconds  Interval  to  GEOS-3 
Revolutions  154,  268,  439,  453  and  467. 


GEOS-3 
Rev.  # 

Raw  R Data 

R.  M.S.  Mean 

Smoothed  R Data 
R.M.S.  Mean 

Residuals  After 
Spline  Fit 
R.M.S.  Mean 

154 

0,2321 

-0.0943 

0.2167 

-0.0943 

0.0833 

0.0000 

268 

1.2643 

1. 1730 

1.2626 

1, 1730 

0.0648 

0.0000 

439 

0,1286 

-0.0049 

0.0933 

-0.0049 

0.0885 

0.0000 

453 

0.4389 

0 . 40  9 6 

0.4333 

0.4096 

0.0695 

0.0000 

467 

0.6512 

0.4345 

0.4772 

0.4345 

0.4431 

0.0000 

Units  are  cm/sec. 


Figure  3.20  is  on  the  same  scale  as  Figures  3.9  to  3,18.  When  we  compare  the 
accelerations  in  revolution  268  (Figure  3,20)  with  accelerations  in  revolutions  154 
and  453  (Figures  3.11  and  3.16),  the  accelerations  R in  revolution  268  appear 
reasonable.  But  there  does  appear  to  be  some  systematic  error  in  R in  revolution 
268  from  Figure  4.19,  when  we  compare  it  with  Figures  3,2  and  3.4.  We  will 
comment  on  this  again  in  Section  4. 1, 

3.3  Radial  Derivative  of  the  Anomalous  Potential  at  GEOS-3  Locations 

The  radial  derivative  ( 9T/3r)c  of  the  anomalous  potential  at  a particular 
location  (<p,  X,  h)  of  the  GEOS-3  satellite  is  obtained  from  the  residual  acceleration 
(R)^  X h by  equation  (1.2)  reproduced  here  specifically: 

(3.32)  = (R)(p,X,h/‘'°®^ 

where  j8  is  the  angle  (Rummel  et  al.,  1976,  page  19)  between  the  radial  direction 
to  GEOS-3  and  the  direction  GEOS-3  to  ATS-6.  Whenever  GEOS-3  happens  to  be 
directly  below  the  ATS-6,  (<p  0°,' X^  266°)  3 T/^  r would  equal  R,  but  at  other 

locations  of  GEOS-3  the  ^T/a  r values  will  be  magnified  sec  j3  times  the  R values. 
We  show  this  for  revolutions  154  and  453  in  Table  3.3  at  one  minute  intervals  out 
of  the  GEOS-3  locations  actually  used  for  the  recovery  of  anomalies  in  Section  4. 
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Table  3.3.  9T/9r  and  R Values  for  GEC)S-3  Revolutions  154  and  453 


Revolution  154 

Revolution  453 

(0° 

■ 

R 

mgals 

sec  j3 

?^T/ar 

mgals 

^0° 

X° 

ii 

mgals 

1 

sec  0 

aT/3r 

mgals 

34 

300 

-0.30 

-1.71 

0.51 

39 

-1.03 

-1.71 

1.75 

31 

297 

-0.26 

-1.57 

0.41 

36 

1.42 

-1.56 

-2.21 

28 

295 

0.21 

-1.45 

-0.30 

33 

287 

2.05 

-1.44 

-2.95 

25 

293 

2.57 

-1.36 

-3.49 

30 

285 

-1.19 

-1.34 

1.59 

21 

291 

-5.45 

-1.28 

6.97 

26 

283 

0.45 

-1.26 

-0.57 

18 

289 

-1.02 

-1.22 

1.24 

23 

281 

1.49 

-1.20 

-1.79 

15 

287 

-2.15 

-1.17 

2.51 

20 

279 

0.22 

-1.14 

-0.25 

12 

285 

-1.12 

0.90 

17 

277 

1.27 

-1.10 

-1.40 

1 

14 

275 

0.15 

-1.07 

-0. 16 

( We  find  that  the  magnification  by  sec  P in  practice  is  not  of  great 

I concern  for  the  area  of  investigation  in  this  report.  However,  when  exceeds, 

say,  60°  ( |sec  ] > 2),  we  may  lessen  the  effect  of  this  magnification  by  smoothing 
the  R data  with  spline  nodes  at  80  seconds  apart  in  that  area,  instead  of  60  seconds. 

. For  this  report,  we  used  spline  nodes  at  uniform  spacing,  but  ti\e  algorithm  in 

Section  3. 1 can  accommodate  uneven  spacing  of  the  nodes  in  view  of  the  scaling  of 
the  data  in  each  interval  according  to  equation  (3. 13).  The  spline  nodes  may  thus 
be  fixed  at  60  seconds  apart  for  data  up  to  60°  around  the  subsatellite  point  of  ATS-6, 
and  80  seconds  apart  after  that.  The  80  seconds  spacing  may  serve  an  adequate 
I balance  between  dampening  of  R because  of  larger  spacing,  and  the  magnification 

of  dT/Sr  because  of  larger  jS. 

Rice  (1969,  Vol.  II,  Chap.  10)  discusses  algorithms  for  the  non-linear 
p problem  of  solving  both  for  the  spline  coefficients,  as  well  as  for  the  spacing  of 

? nodes  by  treating  them  as  variable  nodes.  De  Boor  (1977)  has  a package  of  programs 

i which  may  be  modified  for  this  purpose.  This  approach  was,  however,  not  tried  in 

this  report,  as  the  treatment  of  noisy  R data  by  fixed  spline  nodes  appeared  to  be 
adequate  for  the  current  investigations. 

4.  Recovery  of  5°  Gravity  Anomalies 

With  the  3T/9r  values  from  Section  3,  the  recovery  of  residual  gravity 
anomalies  (Ag)  was  attempted  by  least  squares  collocation  according  to  the  pro- 
cedures described  in  Rummel  et  al.(1976).  Eight  5°  anomalies  were  chosen  which 
were  covered  by  GEOS-3  revolutions  described  in  Section  2.  We  first  present  in 
Section  4. 1 the  results  of  tests  with  GEOS-3  revolutions  in  Section  2. 1,  where 
reliable  values  of  initial  state  vectors  were  available  for  both  satellites.  Next, 
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we  describe  in  Section  4.2  the  tests  with  revolutions  discussed  in  Section  2.3 
where  the  ATS-6  elements  were  not  known  reliably,  but  were  obtained  by  inte- 
gration from  another  epoch.  We  continue  this  inquiry  in  Section  4.3  where  the 
initial  state  vectors  are  obtained  by  integration  for  both  ATS-6  and  GEOS-3 
satellites  instead  of  using  converged  elements.  The  combined  results  using  all 
reliable  observations  are  presented  in  Section  4.4. 

4. 1 Initial  State  Vectors  Available  for  Both  Satellites 

The  GEOS-3  revolutions  available  from  Section  2.1  were  154,  268,  439, 
and  453.  The  location  of  GEOS-3  at  30  seconds  interval  in  these  revolutions  is 
shown  in  Figure  4.1.  The  location  of  the  eight  5°  anomalies  chosen  for  recovery 
are  also  shown.  The  terrestrial  value  of  these  anomalies  Ag , and  their  standard 
deviations,  were  taken  from  Rapp  (1977,  pages  55,  56)  and  the  same  numbering 
system  was  retained.  The  expected  value  of  the  residual  anomalies  were 

obtained  by  subtracting  from  the  anomalies  implied  by  (12,  12)  potential 
coefficients  in  the  GEM  7 set.  (For  details,  see  Rummel  etal.  (1976,  pages  20, 
21).)  The  values  of  the  eight  5°  anomalies  are  given  in  Table  4.1. 


Table  4.1.  particulars  of  Eight  5°  Residual  Anomalies 


Anom. 

# 

(P° 

B 

|H9| 

mgals 

s.  d. 

mgals 

Agpc 

mgals 

Ag' 

mgals 

402 

35 

30 

289 

283 

-34.1 

2.4 

-22.3 

-11.8 

403 

35 

30 

295 

289 

-26.8 

2.5 

-25. 1 

- 1.7 

465 

30 

25 

287 

281 

-15.4 

2.8 

-23.6 

8.2 

466 

30 

25 

293 

287 

-29.8 

2.8 

-30.2 

0.4 

531 

25 

20 

285 

279 

7.3 

2.6 

- 9.9 

17.2 

532 

25 

20 

290 

285 

-44.9 

2.9 

-14.7 

-30.2 

599 

20 

15 

282 

277 

5.8 

1.4 

5.0 

0.8 

600 

20 

15 

287 



282 

- 5.5 

1.3 

- 2.6 

- 2.9 

The  recovered  values  of  anomalies,  was  compared  with  the  expected 
value  of  the  anomalies  by  examining  the  anomaly  discrepancies  € (Ag)  and  the 
correlation  coefficient  p,  which  were  computed  as; 
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e(Agi)  = ^Vi  - E(Agi) 


^ n . n 

(4.2)  p = (7  Agi  E(Ag'i)/n)  / ((7  AVi®/n)^  (V  E(AgoVn)® 


where  n is  the  number  of  anomalies. 

The  3T/9r  values  used  in  the  recovery  of  anomalies  were  30  seconds 
apart,  and  all  available  data  up  to  7? 5 from  the  center  of  the  anomaly  block  was 
considered.  Different  solutions  were  tried,  considering  the  9T/3r  values  to 
have  a standard  deviation  of  0.  5,  1.0,  1.5,  2.0,  2.0  or  5.0  mgals.  (For  details, 
see  Rummel  etal.  (1976,  Sec.  6).)  We  first  report  the  results  with  all  BT/Sr 
values  having  a standard  deviation  of  1.  5 mgals.  The  reasons  for  this  choice 
will  be  ejqilained  later. 

The  statistics  for  the  recovered  anomalies  is  shown  in  Table  4.2 
using  BT/9r  values  in  GEOS-3  revolutions  154,  268,  439  and  453.  Three 
different  solutions  were  tried  with  BT/9r  values  being  computed  for  the  cases 
when  the  spline  nodes  for  fitting  the  R data  were  40  , 60  or  80  seconds  apart. 


Table  4.2.  Statistics  for  Recovered  Anomalies  Using  GEOS-3 
Revolutions  154,  268,  439,  453 


Spline  Node 
Spacing 
(sec. ) 


RMS  Expec. 
Anom.  E(Ag) 
(mgals) 


RMS  Predic. 
Anom.  Ag 
(mgals) 


RMS  Anom. 
Discr.  e(Ag^) 
(mgals) 


Correln. 
coeff.  p 
(eqn.  4.2) 


The  results  showed  a poor  recovery  of  anomalies.  We  recall  from 
discussion  in  the  end  of  Section  3.2  that  there  was  reason  to  suspect  the  R data 
in  revolution  268.  The  recovery  of  anomalies  was  therefore  attempted  from 
only  three  revolutions,  154,  439  and  453,  though  we  notice  from  Figure  4.1 
that  the  data  then  becomes  sparse  over  the  four  eastern  anomalies,  when  we 
omit  revolution  268.  The  statistics  for  the  recovered  anomalies  is  shown  for 


three  solutions  for  dT/3r  values  being  computed  for  spline  node  spacing  for 
fitting  the  R data  as  40,  60  or  80  seconds.  The  standard  deviation  of  5T/9r 
values  was  assumed  to  be  1.5  mgals. 


Table  4.3.  Statistics  for  Recovered  Anomalies  Using  GEOS-3 
Revolutions  154,  439,  453 


RMS  Expec. 
Anom.  E (Ag) 
(mgals) 

RMS  Predic. 
Anom.  A^^ 
(mgals) 

1 

RMS  Anom. 
Discr.  e (Ag”^) 
(mgals) 

Correln. 
coeff.  p 
(eqn.  4.2) 

40 

13.3 

7.7 

10.4 

0.63 

60 

13.3 

7.8 

10.7 

0.60 

80 

13.3 

_ 7.7  j 

11.3 

0.53 

We  first  note  that  in  spite  of  sparse  data  over  four  eastern  anomalies, 
the  anomaly  recovery  is  in  fact  improved  in  Table  4.  3 as  compared  to  Table  4.2, 
which  shows  that  revolution  268  should  be  taken  out  of  the  solution.  We  also 
found  that  though  the  solutions  do  get  worse  with  data  obtained  through  spline 
nodes  at  80  seconds,  and  the  same  was  found  in  separate  tests  for  spline  nodes 
at  30  seconds,  but  there  is  no  noticeable  difference  in  the  solutions  using  spline 
nodes  at  40  or  60  seconds,  the  40  seconds  solution  appearing  slightly  better. 
However,  in  view  of  the  discussion  in  Section  3.2,  we  would  henceforward  use 
the  data  obtained  by  fitting  splines  at  60  seconds  spacing. 

We  now  examine  what  optimum  standard  deviation  should  be  assumed 
for  the  9T/9r  data  used  in  the  recovery  of  anomalies.  During  simulation  studies 
(Rummel  et  al. , Sec  6. 1),  it  was  found  that  a standard  deviation  lower  than  0.5 
mgals  would  cause  instability  in  the  solution  while  a standard  deviation  of  2 mgals 
or  larger  would  dampen  the  solution.  The  results  of  using  different  standard 
deviations  is  shown  in  Table  4.4  for  the  case  covered  in  Table  4.3,  i.e.  using 
data  in  three  revolutions  154,  439  and  453  at  60  seconds  spacing  of  spline  nodes. 

The  solutions  are  unstable  when  the  standard  deviation  of  ?iT/Sr  data 
is  equal  to  or  less  than  1 mgal,  as  evidenced  by  large  RMS  value  of  anomaly 
discrepancy,  c (Ag),  as  well  as  by  the  large  RMS  value  of  recovered  anomalies. 
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Table  4.4.  Statistics  for  Recovered  Anomalies  for  Different 
Assumed  Standard  Deviation  of  J^T/^r.  GEOS-3 
Revolutions  154,  439,  453.  Spacing  of  Spline 
Nodes  60  seconds 


Assumed  s.  d.  of 
3T/9r  data  (mgals) 

) 

RMS  (^) 
mgals 

RMS  ei^') 
mgals 

rn 

1 

1 

0.5 

17.6 

15.1 

0.56 

! 1.0 

10.5 

10.9 

0.60 

1 1.5 

7.8 

10.7 

0.60 

2.0 

6.2 

10.9 

0.59 

3.0 

4.1  • 

11.5 

0.57 

5.0 

2.1 

12.3 

0.54 

(^).  The  RMS  value  17.6  mgals  of  recovered  anomalies  in  Table  4.4  for  the 
standard  deviation  of  0.  5 mgals,  for  example,  is  much  larger  than  the  RMS  value 
of  13.3  mgals  for  the  expected  anomalies.  We  also  note  from  Table  4.4  that  the 
solutions  are  dampened  off  when  the  standard  deviation  is  2 mgals  or  larger,  as 
evidenced  by  smaller  RMS  values  of  the  recovered  anomalies.  The  solution 
appears  to  be  optimum,  when  the  standard  deviation  of  3T/9r  data  is  assumed 
to  be  1.5  mgals.  Hence,  the  statistics  of  recovered  anomalies  will  be  reported 
for  the  standard  deviation  of  1.5  mgals,  as  was  already  done  in  Tables  4.2  to  4.4. 

To  ensure  our  conclusion  (after  Table  4.3)  that  the  data  in  revolution 
268  had  some  systematic  error  and  was  thus  responsible  for  the  poor  recovery 
seen  in  Table  4.2,  we  tried  two  other  solutions  shown  in  Table  4.  5.  In  the  first 
case  (first  row  in  Table  4.5),  we  note  from  Figure  4.1  that  as  revolution  439  was 
not  directly  located  over  the  8 anomalies,  it  could  only  marginally  improve  their 
recovery  (second  row  in  Table  4.3),  as  compared  to  using  revolutions  154  and 
453  only.  Secondly,  if  the  effect  of  revolution  439  is  indeed  marginal,  and  if  the 
data  in  revolution  268  is  erroneous,  the  results  of  using  revolutions  268  and  439 
should  be  much  poorer.  This  was  actually  found  to  be  so,  as  may  be  seen  in  the 
second  row  of  Table  4.5.  Both  solutions  in  Table  4.5  assumed  the  standard 
deviation  of  dT/9r  as  1.5  mgals,  and  the  spacing  of  spline  nodes  was  60  seconds. 
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Table  4.  5.  Statistics  for  Recovered  Anomalies  Using  Only 
2 GEOS-3  Revolutions  out  of  Revolutions  154, 
268,  439,  453 


Revolutions  Used 

RMS  E (Ag) 

RMS  Afe' 

RMS  e (Ag) 

p 1 

(see  Fig.  4. 1) 

mgals 

mgals 

mgals 

154  and  453 

13.3 

8.3 

11.4 

0.53 

268  and  439 

13.3 

9.5 

15.6 

0.10 

We  have  to,  therefore,  take  the  data  in  revolution  268  out  of  the 
solution.  (This  will  be  discussed  again  in  Section  4.3.)  We  would,  of  course, 
prefer  the  solution  using  three  revolutions  154,  439  and  453,  as  shown  in  the 
second  row  of  Table  4.3  over  the  solution  using  two  revolutions  154  and  453,  as 
shown  in  the  first  row  of  Table  4.5.  But,  as  revolution  439  lies  outside  the  eight 
anomalies  to  one  side,  it  contributes  to  the  recovery  only  marginally.  It  is  there- 
fore remarkable  that  we  do  get  a recovery  with  only  two  revolutions  154  and  453. 
The  location  of  data  points  at  30  seconds  interval  used  in  these  two  revolutions 
with  respect  to  eight  anomalies  is  shown  in  Figure  4.2.  In  spite  of  the  data  being 
so  sparse,  the  solution  in  the  first  row  of  Table  4.  5 shows  that  the  procedures 
are  satisfactory.  The  improvement  in  the  solution  shown  in  the  second  row  of 
Table  4.3  with  the  addition  of  revolution  439  is  a pointer  that  the  recovery’  of 
anomalies  would  show  improvement,  as  we  add  to  the  data  set  from  more  revo- 
lutions located  over  the  anomalies. 


4.2  Initial  State  Vector  Obtained  by  Integration  for  ATS-6  Satellite 

As  we  have  already  considered  the  data  in  all  revolutions  available  in 
Section  2. 1,  we  can  only  consider  additional  data  from  the  revolutions  in  Section 
2.3,  where  the  initial  state  vector  was  available  only  for  GEOS-3,  while  that  for 
ATS-6  had  to  be  obtained  by  integration  from  an  earlier  epoch.  VV’e  recall  from  the 
discussion  at  the  end  of  Section  2.3  that  we  had  data  available  for  consideration 
for  three  revolutions,  which  we  called  232  A,  2321  and  2461.  The  initial  ATS-6 
elements  were  obtained  from  integration  from  the  epoch  0“  on  25  April  75  for 
revolutions  2461  and  2321,  and  were  the  adjusted  values  in  revolution  232  A for 
56  minutes  of  range-rate  data.  The  initial  state  vector  for  GEOS-3  was  kept  fixed 
as  the  a priori  value  in  Table  2.  3.  The  RMS  and  the  mean  value  of  the  raw  residual 
range-rate  R data,  the  smoothed  R data  after  fittii^  a spline  with  nodes  60  seconds 
apart,  and  the  residuals  remaining  after  the  spline  fit,  are  shown  in  Table  4.6  for 
the  three  revolutions. 
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Figure  4.2  Location  of  Data  for  Recovery  of  Eight  5°  Gravity  Anomalies 


Table  4.  6.  Least  Squares  Spline  Fit  to  Raw  R 
Data  in  GEOS -3  Revolutions  232  A, 
232 1 and  246  I 


Revolution 

# 



Raw  R Data 

Smoothed  R Data 

Res.  after 

Spline  Fit 

RMS  1 

1 

Mean 

RMS 

Mean 

232  A 

13.3072  j 

-13,1831 

13.3071 

-13. 1831 

0.0643 

0.0000 

232  I 

8. 5816  i 

1 

- 8.5790 

8.5813 

- 8.5790 

0.0643 

0.0000 

246  I 

1, 

2.3917  ' 

- 2.2856 

2.3909 

- 2.2856 

0.0640 

0.0000 

Units  are  cm/sec . 


When  we  compare  the  RMS  and  mean  value  of  the  raw  data  in  Table  4.  6 
with  those  in  the  revolutions  from  Section  2, 1 shown  in  Table  3.2,  we  find  the 
values  for  revolution  246  I to  be  large.  This  may,  however,  be  explained  because 
of  the  uncertainty  in  the  initial  state  vector  of  ATS-6  obtained  by  integration  from 
an  epoch  2 days  earlier.  But  the  much  larger  values  in  revolution  232  I point  to 
some  systematic  error  in  the  raw  R data  in  this  revolution.  This  was  also  pointed 
out  in  the  discussion  after  Table  2.  5.  We  therefore  do  not  expect  a good  recovery 
from  data  in  revolution  232  I. 

The  values  in  Table  4.6  for  the  raw  R data  in  revolution  232  A are 
surprising  when  compared  to  the  values  for  revolution  232  I,  as  the  purpose  of 
letting  the  initial  state  vector  of  ATS-6  be  adjusted  was  to  make  the  residuals  for 
the  full  GEM-7  set  of  potential  coefficients  to  be  small  in  the  least  squares  sense 
(see  items  2 and  3 under  iteration  1 and  7 for  revolution  232  in  Table  2.5).  However, 
in  dispersing  the  large  residuals  (in  iteration  1),  the  ATS-6  initial  state  vector  was 
seriously  distorted  from  the  a priori  integrated  values  (see  Table  2.6)  causing  a 
worsening  of  raw  R data  (residual  to  degree  12  reference  field)  instead  of  improving 
it  in  Table  4.6.  This  is  shown  in  a very  imeresting  manner  by  the  plot  of  acceler- 
ations in  Figures  4.3  and  4.4  for  revolutions  232  I and  232  A respectively.  The 
pattern  of  accelerations  in  revolution  232  A is  the  same  as  in  revolution  232  I,  but 
has  been  severely  biased,  even  running  out  of  the  limits  in  Figure  4.4. 

We  may  recall  here  the  discussion  after  Tables  2.  5 and  2.6.  We  incur 
a serious  risk  in  obtaining  convergence  for  the  initial  state  vector  from  a single 
observation  type  or  for  a limited  time  span.  If  the  observations  are  biased,  the 
converged  initial  state  vector  obtained  on  the  basis  of  fitting  these  observations  in 
the  least  squares  sense  would  be  much  in  error.  We  therefore  do  not  expect  any 
recovery  at  all  from  data  in  revolution  232  A. 
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Data  Interval  10  sec. 


8 mgals 


Figure  4.4  Smoothed  Accelerations  for 
GEOS-3  Revolution  232  A 
with  Spline  Nodes  60  sec. 
apart.  Data  interval  10  sec. 


The  statistics  for  the  recovery  of  6 anomalies  using  data  for  only  one 
revolution  232  A,  232  I and  246  1 at  a time  is  shown  in  Table  4.7.  The  spline 
nodes  were  60  seconds  apart  and  the  standard  deviation  of  3T/3r  was  assumed 
as  1. 5 mgals.  The  north-east  and  the  south-west  anomalies  were  not  considered, 
as  there  would  be  no  data  within  7? 5 from  the  center  of  these  anomaly  blocks. 

This  is  seen  in  Figure  4.5  from  the  location  with  respect  to  eight  anomalies  of 
the  3T/3r  data  points  at  30  seconds  interval  used  in  revolutions  232  and  246. 


Table  4.7.  Statistics  for  6 Recovered  Anomalies  Using  One  GEOS-3 
Revolution  at  a Time  Out  of  Revolutions  232  A,  232  I, 
246  I 


Revolutions  Used 
(see  Figure  4.  5) 

RMS  E(Ag') 
mgals 

RMS(Afe') 

mgals 

RMS  e (Ag) 
mgals 

- - 1 

232  A 

15.4 

26.8 

31.9 

-0.08 

232  I 

15.4 

6.0 

18.0 

-0.29 

246  I 

15.4 

15.2 

14.3 

0.56 

The  data  in  revolution  232  is  unusable,  whether  the  initial  state  vector 
of  ATS-6  was  obtained  by  integration  (232  I)  or  was  then  adjusted  (232  A)  to  fit  the 
56  minutes  of  range-rate  sum  observations.  However,  it  appears  that  the  data  in 
revolution  246  1 is  usable  and  may  help  improve  the  solution  obtained  earlier  with 
j 3 revolutions  154,  439,  and  453  in  Table  4.3.  The  statistics  for  the  recovery  of 

eight  5°  anomalies  are  shown  in  Table  4.8  with  2,  3 and  4 revolutions.  The  spline 
j nodes  in  all  cases  were  60  seconds  apart  and  the  standard  deviation  of  9T/9r  data 

' was  taken  as  1.  5 mgals. 

I 


Table  4.8.  Statistics  for  Recovered  Anomalies  Using  GEOS-3 
Revolutions"  154,  439,  453  and  246  I 


We  find  continuing  improvement  in  the  solutions,  as  data  from  more 
revolutions  is  added,  by  a reduction  in  the  anomaly  discrepancy  e (Ag)  and  increase 
in  the  correlation  coefficient  p (equations  4.1  and  4.2).  But  what  is  more  remark- 
able is  that  the  data  in  revolution  246  I is  usable  inspite  of  the  ATS-6  elements 
having  been  integrated  from  an  epoch  two  days  earlier  (Table  2,  3)  and  thus  causing 
large  residual  range-rate  (Table  4.6).  The  modeling  error  in  K due  to  uncer- 
tainty in  ATS-6  initial  state  vector  was  apparently  'differenced  out'  during  numerical 
differentiation  and  did  not  significantly  bias  the  ii  (and  9T/9r)  values.  However, 
if  this  contention  is  true,  we  should  be  able  to  integrate  from  well  determined  (by 
multiple  observation  types  and  long  time  span)  initial  state  vectors  at  0^“  on 
April  25,  1975  (Table  2.  3)  not  only  for  the  ATS-6  satellite  but  also  for  the  GEOS-3 
satellite.  We  would  then  be  able  to  examine  the  difference  in  integrated  and 
previously  converged  elements  of  GEOS-3  (on  the  lines  of  Table  2.4  for  ATS-6)  and 
check  out  the  values  of  R,  R and  recovered  anomalies  not  only  for  revolutions  232 
and  268  where  the  data  is  suspected  of  having  systematic  errors,  but  also  for 
revolution  246  where  the  range  rate  sum  observations  have  already  yielded  satis- 
factory results. 


4.3  Initial  State  Vectors  Obtained  by  Integration  for  Both  ATS-6  and  GEOS-3 
Satellites 

The  differences  in  the  initial  state  vector  for  GEOS-3  satellite  as 
obtained  from  integration  from  the  epoch  o"  on  April  25,  1975  minus  the  converged 
initial  state  vector  values  in  Tables  2.1  and  2.3  is  shown  in  Table  4.9.  The 
integration  was,  of  course,  carried  out  only  for  the  period  (April  25  to  AprQ  29, 
1975)  during  which  the  observations  were  used  for  determinii^  the  initial  state 
vector  at  O'"  on  April  25,  1975  (last  column  in  Table  2.3).  This  was  discussed 
after  Table  2.4  where  we  concluded  that  the  integrated  elements  could  only  be 
j compatible  during  the  period  on  which  they  were  converged. 

! 


Table  4.9.  Differences  in  Position  and  Velocity  Coordinates  for 
GEOS-3  as  Integrated  From  Epoch  O"  on  April  25, 
1975  Minus  Converged  Elements  for  Other  Epochs 


Epoch 

(Date,  Hr. , min. , sec.) 

X 

(m) 

Y 

(m) 

z 

(m) 

X 

(cm/sec) 

Y 

(cm/sec) 

Z 

(cm/sec) 

1.  26  April  75  07h  00m  00s 
(GEOS-3  revoln.  232) 

-17.2 

50.3 

-22.4 

2.28 

-6.39 

-2.00 

2.  27  April  75  07h  00m  00s 
(GEOS-3  revoln.  246) 

- 7.4 

-1.7 

18.0 

- 0.15 

2.56 

-0.71 

3.  28  April  75  20h  50m  00s 
(GEOS-3  revoln.  268) 

30.3 

573.8 

1 

201.4  ! 

-23.25 

16.44 

40.36 
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The  differences  are  small  for  revolution  246,  though  substantially  what 
may  be  expected  in  converged  elements,  being  10  to  20  m in  position  elements  and 
1 to  2 cm/sec  in  velocity  elements.  The  differences  are  much  larger  for  revolu- 
tion 232,  20  to  50  m in  position  elements  and  2 to  6 cm/sec  in  velocity  elements. 
There  are,  however,  grossly  large  differences  in  the  integrated  and  convei-ged 
elements  of  GEOS-3  for  revolution  268.  We  had  also  noticed  very  large  differences 
for  this  revolution  for  the  integrated  and  converged  ATS-6  elements  (row  8b  in 
Table  2.4).  Some  observations  in  the  2 1/2  hour  time  span  used  for  converging 
revolution  268  in  Table  2. 1 are  likely  to  have  systematic  bias.  Also,  it  is  con- 
ceivable that  when  the  time  span  is  comparatively  shorter  (2  1/2  hours  vs.  5 days), 
a very  different  set  of  initial  state  vectors  for  a pair  of  satellites  may  fit  the  data. 
The  initial  state  vectors  based  on  the  longer  time  span  would  be  more  reliable. 

We  will  refer  to  the  data  in  these  revolutions  as  232  II,  246  II  and 
268  n to  indicate  that  the  initial  state  vectors  of  both  satellifes  have  been  obtained 
by  integration  from  the  epoch  O’’  on  April  25,  1975.  The  particulars  of  R data  in 
these  revolutions  are  given  in  Table  4. 10. 


Table  4. 10.  Least  Square  Spline  Fit  to  Raw  R Data  in  GE06-3 
Revolutions  232  n,  246  II,  268  II 


Revolution 

• 

Raw  R Data 

Smoothed  R Data 

Res.  after 

Spline  Fit 

# 

RMS 

Mean 

RMS 

Mean 

RMS 

— 1 

Mean 

232  II 

5.1172 

-5.1050 

5.1168 

-5. 1050 

0.0643 

0.0000 

246  II 

4.3777 

-4.3495 

4.3772 

-4.3495 

0.0640 

0.0000 

268  II 

1.2855 

-1.1834 

1.2839 

-1.1834 

0.0648 

0.0000 

Units  are  cm/sec . 


The  magnitude  of  R is  reduced  as  compared  to  revolution  232  I while 
slightly  enlarged  as  compared  to  revolution  246  I.  The  latter  should  be  expected 
as  the  integrated  elements  of  GEOS-3  from  a 2 day  earlier  epoch  in  revolution 
246  II  have  greater  uncertainty  and  thus  large  residual  R.  (Compare  with  Table  4.  6). 

The  RMS  values  of  R in  revolution  268  II  are  about  the  same  as  for 
revolution  268  (Table  3. 2),  but  the  mean  values  have  an  opposite  sign.  The  plot 
of  R in  revolution  268  II  is  shown  in  Figure  4.  6.  It  is  on  the  same  scale  and  uses 
the  same  data  as  in  Figure  3. 19  and  is  remarkably  similar,  except  for  its  slope 
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Figure  4. 6 Smoothing  of  Range-Rate  Data  for  GEOS-3 
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reflecting  a negative  bias  while  the  R data  in  revolution  268  in  Figure  3. 19  had  a 
positive  slope.  The  R values  in  revolution  268  n are  shown  in  Figure  4.7  and 
these  are  also  remarkably  similar  to  those  in  Figure  3.20  for  revolution  268, 
except  for  a slight  difference  in  slope.  The  differences  in  Figures  3. 19  and  4.6 
for  R and  in  Figures  3.20  and  4.7  for  R are  only  due  to  different  initial  state 
vectors. 

The  initial  state  vectors  thus  cause  a systematic  long  wavelength  effect 
with  large  variations  in  R,  but  a much  smaller  effect  on  R.  The  plot  of  R for 
revolution  232  II  is  shown  in  Figure  4. 8,  which  is  almost  the  same  as  R in  Figure 
4.3  for  revolution  232  I.  The  difference  is  remarkably  small  considering  that  the 
GEOS-3  elements  are  different  by  20  to  50  m in  position  and  2 to  6 cm/sec  in 
velocity  coordinates.  (We  recall  that  the  ATS-6  elements  were  the  same  in 
revolutions  232  li  and  232  I.) 

The  pattern  of  R appears  to  be  indicative  of  the  anomalous  potential 
sensed  at  GEOS-3  location,  and  is  qualitatively  dependent  on  the  spline  fit  discussed 
in  Section  3.  However,  it  is  quantitatively  affected  by  the  bias  caused  by  incorrect 
initial  state  vectors.  As  the  difference  in  revolutions  246  I and  246  II  in  the 
initial  state  vectors  of  GEOS-3  (there  was  no  difference  in  ATS-6  elements)  is 
more  representative  of  the  type  of  uncertainty  (10  to  20  m in  position  coordinates, 

1 to  2 cm/sec  in  velocity  coordinates),  we  list  in  Table  4. 11  the  actual  differences 
in  3T/Sr  values  for  some  of  the  points  actually  used  in  anomaly  recovery.  The 
points  are  one  minute  apart  in  time. 


Table  4.11.  Variation  in  BT/Br  due  to  Change  in  Initial  State  Vector 
for  GEC)S-3  in  Revolutions  246  I and  246  II  (row  2 in 
Table  4.9) 


<P^ 

BT/Br  in  mgals 

Difference 

Rev.  246  II 

Rev.  246  I 

in  mgals 

(1) 

(2) 

(3) 

(4) 

(5)  = (3)-  (4) 

17 

293 

5.6 

6.4 

-0.8 

20 

291 

3.9 

4.8 

23 

289 

4.8 

5.8 

-1.0 

27 

287 

1.2 

2.3 

-1.1 

30 

284 

-0.3 

1.0 

-1.3 

33 

282 

-1.4 

-1.4 

36 

280 

-0.4 

1.1 

-1.5 

. 
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8 mgals 


We  note  a linear  bias  in  3T/9r  values  due  to  changes  in  initial  state 
vector.  The  differences  in  9T/9r  arc  magnified  bj  about  6 times  in  the  differences 
in  the  recovered  anomalies  using  revolutions  24G  I and  24G  II.  This  maximum 
difference  is  in  the  anomalies  directly  below  the  GEOS-3  locations  and  gets  reduced 
in  neighboring  anomalies.  Accordingly,  the  largest  difference  which  may  be 
expected  in  recovered  anomalies  due  to  uncertainties  in  initial  state  vectors  is 
about  one-half  of  the  standard  deviation  of  the  anomalies,  which  was  about  12  mgals. 

The  statistics  for  the  six  5°  anomalies  recovered  by  data  in  revolutions 
232  II  and  246  II,  one  revolution  at  a time,  is  given  in  Table  4. 12.  These  values 
are  comparable  to  entries  in  Table  4.7.  The  statistics  for  eight  5°  anomalies 
recovered  by  revolutions  268  and  268  II  are  also  given  in  Table  4. 12. 


Table  4.12.  Statistics  for  Recovered  Anomalies  Using  One  GEOS-3 
Revolution  at  a Time  out  of  Revolutions  232  II,  246  II, 
268  II 


Revolutions 

Used 

1 

RMS  E (Ag) 
mgals 

RMS  (4?) 
mgals 

RMS  e(Ag'') 
mgals 

P 

232  II 

15.4 

7.7 

18.2 

-0.14 

246  II 

1 

15«4  I 

11.6 

12.8 

0.58 

268  II 

13.3 

12.0 

14.9 

0.31 

268 

13.3 

9.5 

16.1 

0.03 

We  note  that  the  recovery  is  equally  poor  in  revolutions  232  I and 
232  n.  The  recovery  in  revolution  268  n is  slightly  better  than  in  revolution  268. 
Apparently,  there  are  some  systematic  errors  in  the  range-rate  sum  observ'ations 
in  revolutions  232  and  268,  and  this  data  has  to  be  excluded  from  the  solutions. 

The  recovery  in  revolution  246  II  happens  to  be  slightly  better  than  for 
revolution  246  I.  But  it  is  preferable  anyway,  to  use  the  initial  state  vectors  for 
both  satellites  analogously,  i.e,  integrated  from  the  values  at  epoch  0^  on 
April  25,  1975,  which  W'as  done  for  revolution  246  II. 
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Now  that  we  have  established  that  the  initial  state  vectors  may  be 
obtained  from  integration,  we  would  like  to  see  if  we  may  obtain  R data  for  any 
other  revolution  covering  the  eight  5°  anomalies,  starting  from  the  epochs  at 
O"  on  April  19  and  O**  on  April  25,  1975  (last  two  columns  of  Table  2.  3)  untU 
April  29,  1975.  The  only  range-rate  sum  observations  over  the  anomalies 
during  this  period,  besides  revolutions  already  considered  (154,  175,  232,  246, 
268),  were  in  GEOS-3  descending  revolution  254  (April  27,  1975  22*')  and 
ascending  revolution  260  (April  28,  1975  08*’).  However,  the  location  of 
revolution  254  closely  followed  the  location  of  revolution  453,  and  would  thus 
not  provide  any  additional  information.  The  R data  was,  however,  obtained 
from  revolution  260,  and  we  may  call  it  revolution  260  II,  as  the  initial  state 
vectors  for  both  ATS-6  and  GEOS-3  satellites  were  obtained  by  integration  from 
epoch  0”  on  April  25,  1975.  The  solution  with  the  additional  data  in  revolution 
260  II  will  be  reported  in  Section  4. 4. 


4.4  Recovery  of  5°  Anomalies  — Combined  Solution 

We  have  finally  the  R data  available  for  revolutions  154,  439,  453, 
246  n and  260  II.  The  location  of  GEOS-3  satellite  at  30  seconds  time  interval 
in  these  revolutions  used  for  the  recovery  of  eight  5°  anomalies  is  shown  in 
Figure  4.9.  The  3T/Sr  values  in  the  locations  where  the  ascending  and  de- 
scending GEOS-3  revolutions  cross  each  other  are  given  in  Table  4. 13. 

Though  the  locations  (o,  X , h)  of  the  crossover  points  were  not 
corresponding  exactly,  but  the  data  could  perhaps  be  adjusted  in  a linear 
manner  in  all  ascending  and  descending  revolutions  so  that  the  discordance  at 
the  crossover  points  is  a minimum  in  the  least  squares  sense.  The  justification 
for  the  linear  adjustment  of  data  in  each  arc  lies  in  the  fact  that  the  errors  in  the 
initial  state  vectors  cause  a linear  bias  in  the  3T/Br  values  (see  Table  4.11). 

The  adjustment  was  not  done  in  the  present  study,  as  we  should  first  perhaps  use 
R data  from  at  least  2 more  revolutions  — one  ascending  revolution  in  a location 
corresponding  to  revolution  232  (Figure  4.5)  and  one  descending  revolution  in  a 
location  corresponding  to  revolution  268  (Figure  4. 1).  Preferably,  we  should 
consider  data  from  another  two  ascending  revolutions,  one  each  in  the  north-east 
and  south-west  comers,  so  that  each  anomaly  is  covered  by  at  least  one  ascending 
and  one  descending  revolution. 

^ We  now  give  in  Table  4. 14  the  values  of  the  anomaly  discrepancies, 

for  each  of  the  eight  5°  anomalies  as  they  change  from  a solution  using 
data  for  only  2 revolutions  (154  and  453),  and  as  data  from  revolutions  439,  246  II 
and  260  II  is  added  for  one  revolution  at  a time.  The  RMS  value  of  the  predicted 
anomalies  (RMS  Ag)  and  the  anomaly  discrepancies  (RMS  (c  Ag)),  and  the  corre- 
lation coefficient  (p)  is  also  given  for  ease  of  reference.  We  recall  that  the  RMS 
value  of  the  expected  anomalies  was  13.3  mgals. 
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Figure  4.9  Location  of  Data  for  Recovery  of  Eight  5°  Gravity  Anomalies 


Table  4.13.  Comparative  3T/9r  Values  at  Crossover  Points  of  Ascending  GEOS-3 
Revolutions  246  II  and  260  II  with  Descending  GEOS-3  Revolutions 


1 


Table  4. 14.  Improvement  in  Anomaly  Recovery  with  R Data 
in  Additional  GEOS-3  Revolutions 


Anom.  # 

Latitud. 

Longitud. 

m mgals  with  R Data  in  Revolns. 

(see  T. 4.1) 

Extent 

Extent 

154  + 453 

+ 439 

+ 246  II 

+ 260  II 

402 

30-35 

283-289 

21.2 

18.2 

18.6 

14.7 

403 

30-35 

289-295 

10.8 

10.5 

10.4 

4.3 

465 

25-30 

281-287 

-9.6 

-7.5 

-9.7 

-9.3 

466 

25-30 

287-293 

5.2 

5.2 

0.3 

-2.3 

531 

20-25 

279-285 

-10.2 

-11.6 

-14.3 

-14.3  j 

532 

20-25 

285-290 

14.7 

14.7 

0.8 

0.6  1 

599 

15-20 

277-282 

4.1 

4.7 

4.7 

4.7  ' 

600 

15-20 

1_ 

282-287 

-3.4 

-3.7 

-3.7  1 

RMS  Anom.  Discr.  ( c ( Ag))  mgals 

11.4 

10.7 

9.9 

8.4  ! 

RMS  Pred.  Anom. 

(Ag)  mgals 

8.3 

7.8 

11.6 

11.0 

Correln 

. Coeff.  (p) 

0.53 

0.60 

0.69 

0.78  ' 

Mean  of  std.  dev.  of  pred.  anom. 
mgals 

12.8 

12.7 

12.4 



12. 1 

We  note  from  Table  4. 14  and  Figure  4.9  that  the  data  in  any  revolution 
primarily  affects  the  anomaly  directly  below  it.  Hence,  as  already  mentioned,  we 
should  have  R data  from  at  least  two  (preferably  four)  additional  revolutions  to 
further  improve  the  anomaly  recovery  in  the  last  column  of  Table  4. 14.  We  should 
also  attempt  to  adjust  the  9T/9r  data  linearly  in  all  revolutions  for  a minimum 
variance  of  non-agreement  over  the  cross-over  points  before  the  final  solution  to 
cater  for  systematic  biases  in  9T/9r  data  due  to  inaccurate  initial  state  vectors. 


Finally,  the  picdicted  anomalies  with  the  five  revolutions  in  Figure  4.9 
used  in  the  present  study  (without  any  cross-over  constraints),  along  with  standard 
deviations,  are  shown  in  Table  4.15.  The  number  of  9T/9r  points  within  3? 5 
{(i)  = 3? 5)  of  the  center  of  each  anomaly  block  are  also  given,  as  these  primarily 
affect  the  anomaly  recovery.  All  9T/9r  values  lying  within  7?5  {ip=  7°.  5)  of  the 
center  of  each  anomaly  block  were,  however,  used;  their  number  is  also  given. 

The  time  interval  of  3T/9r  values  was  30  seconds  and  they  were  assumed  to  have 
a standard  deviation  of  1.5  mgals  to  ensure  stability  of  the  solution  (see  Table  4.4 
and  relevant  discussion).  The  splines  to  filter  and  smooth  the  raw  R data  had 
their  nodes  60  seconds  apart,  and  used  R data  at  10  seconds  interval. 
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Table  4.15.  Anomaly  Recovery  of  Eight  5°  Anomalies  Using  R Data 
From  Five  GEOS-S  Revolutions  154,  439,  453,  246  II, 
260 II  (see  Fig.  4.9) 


Anomaly 

Latitud. 

Longitud. 

Expec. 

# of  Data 

Predic. 

1 

Stnd. 

Anom. 

# 

Extent 

Extent 

Anom. 

Points 

Anom. 

Dev. 

Discr. 

mgals 

l/)=3^5  ([)=7?5 

mgals 

mgals 

mgals 

402 

30-35 

283-289 

-11.8 

9 

30 

2.9 

11.3 

14.7 

30-35 

289-295 

- 1.7 

3 

25 

2.6 

12.4 

4.3 

465 

25-30 

281-287 

8.2 

27 

-1.1 

11.3 

- 9.3 

25-30 

287-293 

0.4 

28 

-1.9 

11.6 

- 2.3 

531 

279-285 

17.2 

18 

2.9 

11.9 

-14.3 

285-290 

-30.2 

25 

-29.6 

12.4 

0.6 

599 

277-282 

0.8 

11 

5.5 

12.7 

4.7 

282-287 

- 2.9 

16 

-6.6 

13.5 

- 3.7 

RMS  Value  (mgals) 

13.3 

11.0 

8.4 

Correlation  coefficient  between  predicted  and  expected  anomalies  = 0.78 


A solution  was  also  tried  with  the  above  data  but  with  standard  deviation 
of  9T/9r  data  as  1.0  mgal,  but  the  solution  was  worse.  The  RMS  value  of  pre- 
dicted anomalies  was  15.7  mgals,  and  the  RMS  value  of  anomaly  discrepancies  was 
10. 7 mgals,  with  correlation  coefficient  as  0.74  mgals.  Additional  data  from 
revolution  268  II  (row  3 of  Table  4. 12)  also  made  the  solution  worse.  The  solution 
in  Table  4. 15  is  therefore  the  optimum  one  with  the  data  from  five  revolutions 
shown  in  Figure  4.9. 

5.  Summary  and  Conclusions 

There  are  two  main  considerations  in  the  use  of  real  range-rate  sum 
observations,  from  which  we  subtract  the  computed  value  of  the  range-rate  sum 
based  on  a certain  model  to  obtain  line  of  sight  residual  range-rate  between  the 
'high-low'  satellite  pair.  First  are  the  modeling  errors  most  significantly 
contributed  by  the  assumed  initial  state  vectors  of  the  satellites,  and  second  are 
the  random  observational  errors  of  the  raw  data,  after  it  has  been  corrected  for 
physical  and  instrumental  effects  like  ionospheric  and  tropospheric  refraction, 
transponder  delays,  etc.  The  observations  may  sometimes  also  exhibit  system 
biases,  which  may  be  accounted  for  if  'small'  and  linear,  but  which  otherwise 
may  make  the  observations  unusable. 
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The  initial  state  vectors  are  'determined'  on  the  basis  that  if  the  earth's 
gravitational  field  at  the  satellite  altitude  may  be  described  with  neglibible  error  by 
a full  set  (e.g.  GEM-7)  of  potential  coefficients  (and  if  we  consider  other  force 
fields  like  solar  radiation  pressure,  etc.),  then  the  computed  values  should  fit  the 
observations  in  the  least  squares  sense  with  the  only  parameters  being  the  initial 
state  vectors.  However,  in  practice,  this  depends  on  the  observation  type  and 
its  reliability  and  the  time  span  of  observations.  Range-rate  sum  observations  are 
not  sensitive  enough  to  the  determination  of  the  initial  state  vectors,  and  if  they 
are  used  on  a 'single  pass’  basis  — the  time  duration  of  the  'visibility'  from  ATS-6 
in  one  revolution  of  GEOS-3— , the  initial  state  vectors  even  if  'converged'  may 
be  seriously  biased. 

As  we  are  interested  in  the  range-rate  sum  observation  type  for  the 
recovery  of  gravity  anomalies,  we  require  a two  stage  process,  the  first  stage 
being  the  determination  of  unbiased  initial  state  vectors  in  which  we  employ  other 
observation  types  particularly  range  observations,  including  laser  ranges,  in 
addition  to  range-rate  sum  observations.  It  is  usually  possible  with  this  multiple 
observation  type  to  get  converged  initial  state  vectors  on  a single  pass  basis. 
However,  if  the  observations  include  sufficient  'paired'  observations  like  range- 
rate  sum  or  range  sum,  one  may  obtain  a 'wrong'  pair  of  initial  state  vectors 
which  may  or  may  not  introduce  a bias  in  range  rate  residuals.  It  appears  to  be  a 
much  better  procedure  to  use  converged  initial  state  vectors  for  a much  longer 
time  span  like  4 to  6 days  (Table  2.  3).  Firstly,  any  biased  observations  in  one 
or  two  revolutions  would  not  significantly  bias  the  initial  state  vectors.  This  will 
result  in  greater  reliability  in  terms  of  the  standard  deviations  of  the  initial  state 
vectors,  which  otherwise  on  a single  pass  basis  even  with  multiple  observations 
may  reach  very  large  values  of  a kilometer  in  position  and  several  meters /sec 
in  velocity  (Table  2. 1),  making  it  impossible  to  judge  their  reliability.  Secondly, 
on  a long  time  span  of  4 to  6 days,  it  would  become  possible  to  use  even  range- 
rate  sum  observations  from  multiple  passes  to  obtain  convergence,  which  even  if 
it  requires  several  iterations,  has  the  practical  advantage  of  not  requiring  access 
to  other  observation  types. 

We  find  that  there  is  no  cause  for  concern  for  obtaining  initial  state 
vectors  for  other  epochs  by  integration,  but  these  epochs  must  be  within  the  time 
span  used  for  obtaining  the  converged  elements  (Table  2.  4).  Any  modeling  error 
due  to  incorrect  initial  state  vector  appears  to  show  up  as  a linear  bias  in  the 
residual  range-rate  R and  the  residual  accelerations  R.  This  may  have  a large 
effect  on  R,  but  gets  differenced  out  while  computing  R,  with  very  small  linear 
bias  on  accelerations  (Table  4. 11)  for  the  usual  uncertainty  of  10  to  20  meters  in 
position  and  1 to  2 cm/sec  in  velocity  of  GEOS-3  initial  state  vector. 

It  is  not  possible  to  filter  any  residual  modeling  errors  and  the  obser- 
vational errors  by  letting  the  initial  state  vectors  take  the  'slack'  by  letting  them 
vary  within  the  standard  deviations  (computed  for  them  in  the  first  stage),  while 
computing  residual  range-rate  R.  The  R function  can,  however,  be  filtered  of 
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observational  'noise'  by  approximating  it  in  the  least  squares  sense  by  a piecewise 
continuous  cubic  spline,  with  continuous  first  and  second  derivatives.  It  is 
remarkable  that  bases  exist  for  this  spline  in  terms  of  two  simple  cubic  poly- 
nomials, such  that  a linear  combination  with  only  four  non-zero  coefficients  is 
required  to  represent  the  spline  between  consecutive  nodes,  where  the  adjacent 
piecewise  continuous  cubics  meet.  This  results  in  a simple  computationally 
stable  algorithm  described  in  Section  3.1. 

A unique  approximating  spline  gets  defined  once  the  spline  nodes  are 
selected.  As  the  spacing  bet\\’een  nodes  is  increased,  the  slopes  of  the  spline,  ii 
in  our  case,  get  smoother.  The  optimum  representation  of  R,  in  the  sense  that 
the  variation  'a  nei^boring  points  at  10  to  30  seconds  interval  is  what  may  be 
expected  at  the  satellite  altitude,  occurs  with  spline  nodes  60  seconds  apart  when 
the  raw  R data  is  10  seconds  apart.  A smaller  spacing  of  spline  nodes  causes  R 
to  vary  too  sharply,  while  a larger  spacing  causes  R to  be  dampened.  The  obser- 
vational noise  in  the  revolutions  used  in  this  study  was  found  to  varj"  from  0.06  to 
0.09  cm/sec,  while  in  one  revolution,  which  was  not  used,  it  was  as  lai^e  as 
0.44  cm/sec  (Table  3.2). 

The  spline  should  be  fitted  to  include  raw  R data  2 to  4 minutes  at  each 
end  beyond  the  data  span  needed  for  recovery  of  anomalies,  so  that  the  spurious 
accelerations  at  the  ends  of  data  span  may  be  discarded.  The  3T/9r  values  are 
magnified  from  R values  by  secant  times  the  angle  at  GEOS-3  between  the  radial 
directions  and  line  of  sight  ATS-6/GEOS-3.  WTien  this  angle  exceeds  60°,  it  may 
be  advisable  to  make  the  spline  nodes  80  seconds  apart  to  reduce  sharp  variations 
in  ST/Sr.  An  alternative  approach  would  be  to  use  an  algorithm  which  chooses  an 
optimum  spacing  for  the  nodes  — variable  nodes  instead  of  fixed  nodes.  Tests 
need  to  be  carried  out  to  investigate  if  in  our  case  of  uniform  interval  data  points 
representing  a smooth  continuous  R function,  the  non-linear  algorithm  of  variable 
nodes  will  give  us  an  advantage,  particularly  for  the  case  when  the  ST/9r  values 
are  magnified  more  than  2 times  theR  values. 

The  choice  of  low  degree  reference  field  represented  by  (12,  12) 
potential  coefficients  out  of  GEM-7  appears  to  be  very  suitable  for  computing  R. 

As  these  coefficients  are  well  determined,  the  residual  range-rate  R does  not  have 
any  significant  bias  due  to  errors  in  potential  coefficients.  Secondly,  a degree  12 
field  is  adequate  at  GEOS-3  altitude  to  make  the  covariance  function  of  ST /Sr 
neglibibly  small  beyond  7? 5 thereby  allowing  a relatively  few  number  of  points  at 
30  seconds  data  interval  (Table  4.15)  to  recover  the  anomalies. 

The  improvement  in  the  recovery  of  anomalies  as  data  from  additional 
revolutions  is  added  (Table  4. 14)  indicates  that  we  should  have  data  from  at  least 
one  ascending  and  one  descending  revolution  over  each  anomaly.  As  we  have 
found  that  initial  state  vectors  may  be  obtained  by  integration  from  converged 
elements  at  epochs  4 to  6 days  apart,  it  is  comparatively  easy  to  reduce  data  from 
more  revolutions  to  achieve  the  required  density. 
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The  ST/Sr  data  should,  however,  be  first  examined  for  each  revolution 
sepai'ately  by  predicting  a few  anomalies  in  common  to  determine  if  there  is  a 
system  bias  in  the  range-rate  sum  observations  in  any  revolution.  We  found  two 
such  revolutions  (232  and  268)  in  the  present  investigations. 

The  selected  revolutions  would  still  have  a small  linear  bias  in 
9T/9r  values  due  to  modeling  errors  caused  by  residual  errors  in  the  initial  state 
vectors.  We  need  to  determine  two  unknowns  for  a zero  offset  and  a slope  coeff- 
icient in  each  revolution  to  adjust  the  BT/9r  values  to  have  a minimum  variance 
discrepancy  at  crossover  points  between  ascendit^  and  descending  revolutions. 

The  problem  becomes  'over-determined'  for  more  than  four  ascending  and  four 
descending  revolutions.  As  we  did  not  have  a sufficient  number  of  revolutions 
reduced  for  the  present  study,  the  adjustment  for  the  cross-over  constraints  was 
not  attempted. 

With  the  unadjusted  values  of  5T/9r  values,  we  found  that  the  solutions 
for  the  predicted  anomalies  were  unstable  for  assumed  standard  deviations  of 
ST/3r  less  than  1.  5 mgals.  It  is  likely  that  this  may  be  reduced  with  adjusted 
9T/Sr  values.  This  would  then  result  in  better  agreement  of  predicted  anomalies 
with  the  expected  values  in  terms  of  RMS  values,  and  also  reduced  standard 
deviations  of  the  predicted  anomalies. 

The  final  solution  for  eight  5°  anomalies,  using  range  rate  sum 
observations  in  five  revolutions  shown  in  Figure  4.9,  is  given  at  Table  4.15. 
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